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PREFACE 

The  National  Environmental  Satellite  Service  (NESS)  operates 
a  central  facility  for  determining  accurate  earth-location  in- 
formation for  the  images  derived  from  its  Geostationary  Opera- 
tional Environmental  Satellite  (GOES) .  The  facility  is  called 
the  VISSR  Image  Registration  and  Gridding  System  (VIRGS) . 

The  main  purpose  of  this  document  is  to  set  forth  a  mathe- 
matical formulation  of  the  entire  earth-location  process  for  a 
spin-stabilized  geosynchronous  satellite.   The  implementation 
of  this  process  on  the  VIRGS  is  described.   Preceding  the  main 
topics  is  a  discussion  of  the  relevant  aspects  of  the  GOES  sys- 
tem.  An  error  analysis  and  some  actual  results  are  presented 
in  the  final  section. 

The  original  objective  of  the  VIRGS  was  that  it  provide  a 
stand-alone  capability  to  operationally  produce  highly  accurate 
earth-location  parameters.   The  "stand-alone"  feature  is  impor- 
tant because  it  eliminates  the  expense  of  operating  ranging 
stations  and  performing  separate  data  processing  tasks.   "High 
accuracy"  is  as  important;  without  it  the  users  of  GOES  imagery 
are  faced  with  the  undesirable  alternatives  of  either  using 
poorly  located  data  or  determining  accurate  earth-location 
parameters  on  their  own. 

The  objectives  of  the  VIRGS  have  not  yet  been  fully  real- 
ized.  This  is  due  in  large  to  the  lag-time  between  the  devel- 
opment of  a  new  technique  and  its  full  operational  implementa- 
tion.  A  technique  that  uses  the  positions  of  stars  in  the 
VISSR  image  to  determine  some  of  the  earth-location  parameters 
has  not  yet  been  implemented.   The  star  position  technique 
streamlines  the  earth-location  process  and  provides  the  needed 
accuracy  and  stand-alone  features.   Also  remaining  to  be  fully 
implemented  is  the  closed-loop  capability  of  the  VIRGS--a  capa- 
bility for  complete  quality  control  of  the  earth-location 
process. 

The  authors'  intent  is  to  show  that  the  objectives  of  the 
VIRGS  can  be  achieved;  and  it  is  their  hope  that  the  demonstra- 
tion given  herein  will  aid  in  that  achievement. 

Some  acknowledgements  are  deserved.   Hank  Schmidt  began 
casually  reading  an  early  version  of  the  document  and  ended  up 
making  significant  improvements  to  the  presentation.   Jim 
Ellickson  and  Matt  Jurotich  reviewed  an  early  version.   Ron 
Gird  is  thanked  for  his  encouragement  and  energetic  participa- 
tion in  the  VIRGS  implementation.   Vince  Oliver  and  others, 
inside  and  outside  of  NESS,  who  as  users  actively  supported 
improvements  to  the  operational  earth-location  of  VISSR  images, 
were  instrumental  in  bringing  about  VIRGS.   The  McIDAS  team  at 
the  University  of  Wisconsin  (SSEC)  planted  the  seed  and  then 
provided  their   know-how  in  delivering  the  basic   VIRGS.    Eric 
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Smith,  J.T.  Young  and  others  were  responsible  for  the  develop- 
ment and  demonstration  of  the  early  navigation  software  on 
McIDAS. 

VIRGS  is  operated  and  maintained  daily,  in  accordance  with 
procedures  and  guidelines  passed  through  management,  by  a  group 
consisting  of  satellite  meteorologists  (oddly  enough) ,  software 
analysts,  computer  operators,  and  electronic  technicians. 
Members  of  this  group  are  recognized  for  their  daily  efforts. 
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ABSTRACT.   Aspects  of  the  GOES  system  rele- 
vant to  earth-location  are  reviewed.   The 
essential  features  of  an  earth-location  capa- 
bility that  satisfies  the  needs  of  users  of 
the  data  are  specified.   A  mathematical  for- 
mulation of  the  geometry  and  orbit/attitude 
determinations,  including  the  actual  algo- 
rithms, is  presented.   The  implementation  of 
these  methods  on  VIRGS  is  described.   Final- 
ly, the  rationale  for  a  centralized,  once- 
for-all  users  earth-location  service  and  an 
error  budget  with  some  actual  results  are 
presented. 


I.   INTRODUCTION 

The  Geostationary  Operational  Environmental  Satellites  (GOES) 
are  examples  of  spin-stabilized  geosynchronous  satellites.   The 
GOES  provide  continuous  viewing  of  the  weather  features  over 
the  Americas  and  their  adjacent  waters.   The  images,  or  views, 
are  obtained  every  half-hour,  and  sometimes  more  frequently, 
from  each  of  two  satellites.   The  imaging  instrument  of  the 
GOES  is  called  a  Visible  and  Infra-red  Spin-Scan  Radiometer 
(VISSR) .   About  a  dozen  users  receive  VISSR  data  directly  from 
the  GOES  and  over  a  hundred  others  are  served  by  a  Central  Data 
Distribution  Facility. 

The  GOES  program  began  in  1974  and  is  expected  to  remain  in 
operation  throughout  the  1980's.   During  the  next  ten  years  the 
number  of  uses  and  users  of  GOES  data  is  expected  to  increase. 

In  order  for  the  VISSR  imagery  to  be  useful,  it  must  be 
earth-located.   "Earth-location"  means  different  things  to 
different  users.   It  can  mean  (1)  the  ability  to  determine 

1 


precisely  the  earth-location  (latitude  and  longitude)  of  any 
feature  in  an  image;  (2)  the  presence  of  fixed  surface  features 
of  earth  in  time-lapsed  (animated)  sequences  of  images;  or  (3) 
having  geographic  and  political  boundaries  accurately  laid  over 
the  images. 

Today  many  users  of  VISSR  images  use  manual  techniques  to 
shift  and  align  images.   These  methods  are  cumbersome,  time 
consuming,  and  inaccurate.   They  provide,  at  best,  partial 
solutions  to  the  earth-location  problem.   A  general  solution 
--one  that  enables  a  user  to  map  any  picture  element  (VISSR 
f ield-of-view)  into  precise  earth  coordinates--must  be  based  on 
a  complete  mathematical  formulation  of  the  problem.   There  are 
two  reasons  for  this:   (1)  the  basic  geometry  that  describes  a 
mapping  or  transformation  from  satellite/instrument  (GOES/ 
VISSR)  coordinates  to  earth  coordinates  is  complex;  and  (2)  the 
dynamics  of  the  satellite's  orbit  and  attitude  must  be  modelled 

With  the  information  needed  to  perform  accurate  earth- 
location  inserted  into  the  VISSR  data  stream,  all  direct  VISSR 
users  would  benefit  by  not  having  to  determine  it  on  their  own. 


II.   GENERAL  REVIEW 
A.   The  Satellite 

A  satellite  in  an  ideal  geosynchronous  orbit  remains  motion- 
less with  respect  to  an  earth-based  observer.   For  this  to 
occur,  the  satellite  must  be  in  the  equatorial  plane  and  the 
centrifugal  force  of  its  orbital  motion  must,  at  all  times, 
equal  the  gravitational  force  between  it  and  the  earth.   How- 
ever, there  are  perturbating  forces  that  make  it  highly  imprac- 
tical to  maintain  such  an  ideal  orbit. 

The  GOES  orbit  is  not  perfectly  circular  and  it  is  slightly 
inclined  to  the  equatorial  plane.   Furthermore,  the  perturba- 
tive  forces,  caused  by  anomolies  in  the  earth's  gravitational 
force  field,  the  gravitational  affect  of  the  moon  and  the  sun, 
and  the  pressure  of  the  solar  wind  continually  vary  the  orbit. 
Departure  from  an  ideal  geosynchronous  orbit  causes  the  satel- 
lite subpoint  to  trace  out  a  "distorted  figure  eight"  during 
its  one-day  orbital  period. 

The  VISSR  is  rigidly  mounted  in  an  upright  position  to  an 
axis  about  which  the  satellite  spins  at  100  rpm.   Ideally  the 
attitude  of  the  spin  axis  would  be  normal  to  the  orbit  plane. 
An  ideal  orbit  and  spin  axis  attitude  would  produce  the  same 
"perspective"  in  all  VISSR  images.   In  practice  external 
torques  from  solar  radiation  and  impulses  from  VISSR  scan 
stepping  cause  the  satellite  spin  axis  to  wobble  around  a 
vector  normal  to  the  orbit  plane  in  a  manner  governed  by  rigid 
body  dynamics.   The  practical  manifestation  is  a  nearly  linear 
precession  of  the  spin  axis  over  a  period  of  several  days 
accompanied  by  a  small  nutation  with  a  cycle  period  of  several 
seconds. 

A  set  of  orbit  parameters  determines  the  satellite's  position 
and  a  set  of  attitude  parameters  determines  the  orientation  of 
the  satellite's  spin  axis  relative  to  the  earth,  as  functions 
of  time.   The  satellite's  position  and  orientation  are  control- 
led so  as  to  stay  within  certain  bounds  and  hence  avoid  exces- 
sive apparent  earth  motion  in  the  VISSR  images.   The  method  of 
control  is  to  fire  the  on-board  rockets  in  a  manner  that  appro- 
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priately  alters  the  orbit  and  attitude.   These  are  called 
"maneuvers"  and  such  maneuvers  are  vital  in  maintaining  the 
quality  of  the  GOES  service.   When  these  maneuvers  occur,  sig- 
nificant discontinuities  in  the  orbit  and  attitude  result.   The 
maneuvers  are  usually  scheduled  several  days  in  advance  and  can 
be  made  as  frequently  as  once  a  week.   The  orbit  and  attitude 
discontinuities  caused  by  satellite  manuevers  interrupt  the 
process  of  predicting  the  orbit  and  attitude  state  with  a 
mathematical  model.   Thus  the  earth  location  capability  is  also 
interrupted.   A  navigational  system,  such  as  VIRGS,  requires 
accurate  orbit  and  attitude  states  within  several  hours  fol- 
lowing a  maneuver. 

B.   The  VISSR 

The  imaging  instrument  on  GOES  is  the  VISSR.   The  VISSR  is 
rigidly  mounted  to  the  satellite's  body.   An  image  is  formed  by 
the  VISSR  Field  of  View  (FOV)  scanning  the  earth,  west  to  east, 
as  the  satellite  spins.   After  each  spin,  a  scan  mirror  is 
moved  (stepped  down  one  notch)  so  that  the  next  scan  sweep  is 
slighly  south  of  the  previous  sweep. 

Electronic  sampling  generates  equal  angle  (or  time)  spacing 
between  successive  image  elements  on  each  scan  line. 
Throughout  each  scan  line  the  scan  mirror  is  at  a  fixed 
position.   The  scan  mirror  stepping  increment  is  192 ^radians 
and  the  sampling  interval  is  84/\radians  for  the  infrared 
channel(s)  and  21^  radians  for  the  visible  channel. 

There  is  a  single  infra-red  detector.   It  has  a  192  by  84 
//cradian  field  of  view  and  provides  7  by  3  km  subpoint 
resolution.   There  are  8  vertically  aligned  visible  detectors 
which  are  sampled  simulataneously  giving  a  21/4iradian  square 
field  of  view,  for  0.8  km  resolution  at  the  subpoint. 

The  normal  operation  of  the  VISSR  is  to  acquire  a  full  earth 
disc  image  frame  (1821  successive  scans)  once  every  30 
minutes.   Sometimes  the  VISSR  is  commanded  to  acquire  less  than 
a  full  earth  disc  image  by  scanning  over  any  interval  of  scan 
positions  within  the  range  of  1  to  1821.   When  this  is  done, 
images  can  be  obtained  for  example  during  3,  7,  and  15  minute 
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intervals,  depending  on  the  size  of  the  images.   The  earth- 
location  capability  described  herein  is  compatible  with  these 
various  imaginq  modes. 

The  mechanical  alignment  of  the  VISSR  to  the  satellite  spin 
axis  and  the  sun  pulse  detector  (see  next  paragraph)  is  imper- 
fect.  This  creates  significant  values  for  three  misalignment 
angles  that  are  analogous  to  the  roll,  yaw,  and  pitch  angles  of 
a  spacecraft.   The  VISSR' s  mechanical  frame  could  be  tilted 
forward  (pitch) ,  tilted  to  the  side  (yaw)  as  well  as  being 
rotated  around  the  vertical  (roll).   The  determination  of  these 
misalignment  angles,  which  are  constant  over  periods  of  weeks, 
is  a  vital  step  in  achieving  accurate  earth  location. 

A  sun  sensor  and  an  earth-edge  detector  are  used  to  initiate 
the  sampling  along  each  scan  line.   During  each  spin  of  the 
satellite  the  delay  between  sensing  the  sun  and  detecting  the 
earth's  west  edge  is  measured.   The  delay  that  is  measured  on 
one  spin  n  is  used  to  initiate  sampling  on  spin  n+1.   This 
process  merely  ensures  that  the  earth  view  is  contained  in  the 
data  acquisition  interval.   The  task  of  precisely  referencing 
the  west  border  of  the  image  frame  on  successive  scan  lines  is 
performed  on  the  ground  in  real-time  by  the  Synchronous  Data 
Buffer. 

C.   Functions  of  the  Synchronous  Data  Buffer  (SDB) 

The  VISSR  data  are  transmitted  in  real-time  at  28  Mbps*  to  a 
Command  and  Data  Acquisition  (CDA)  facility  on  the  ground.   At 
the  CDA  an  SDB  is  used  to  "stretch"  the  data  stream  for  re- 
transmission through  GOES  at  1.74  Mbps.   The  stretching  opera- 
tion uses  the  340°  portion  of  each  satellite  revolution,  when 
the  VISSR  is  not  pointing  at  the  earth,  to  rebroadcast  the 
VISSR  data  at  1/16  the  raw  data  rate.   Users  of  the  data  are 


*Megabits  per  second  (Mbps) 


therefore  interested  in  the  form  of  the  SDB's  VISSR  output 
rather  than  the  satellite's  raw  output. 

The  SDB  accomplishes  other  critical  tasks  in  the  process  of 
stretching  the  data  stream.   From  the  point  of  view  of  earth- 
locating  the  data,  one  of  the  SDB's  most  critical  functions  is 
to  center  each  VISSR  scan  line  on  the  earth.   The  SDB,  while 
maintaining  synchronization  with  the  satellite's  spin  rate, 
precisely  times  the  start  (west  border)  of  each  stretched  scan 
line.   To  do  this  the  time  of  the  sun  pulse  and  a  "beta"  angle 
are  used.   The  "beta"  angle  is  the  angle  subtended  at  the 
satellite  by  lines  projected  from  the  sun  and  earth  to  the 
satellite,  in  the  satellite's  spin  plane.   The  beta  angle  is 
computed  by  a  sun  position  model  using  the  satellite's  orbit 
and  attitude  and  the  VISSR  misalignment  parameters.   This 
computation  is  readily  available  (as  a  by-product)  from  an 
earth-location  capability  such  as  VIRGS.   The  value  of  beta 
used  on  each  scan  line  references  each  sample  of  that  scan  line 
with  respect  to  the  sun's  position.   The  well-defined  relation- 
ship between  beta  and  every  sample  of  the  scan  line  provides  an 
absolute  reference.   Hence,  even  if  the  "beta"  provided  to  the 
SDB  is  inaccurate,  the  user  of  the  data  can  still  reference  the 
imagery  data  to  the   sun's  position.   The  "beta"  used  by  the 
SDB  is  documented  in  the  stretched  VISSR  data  stream  (SDB 
output) . 

The  image  frame  as  seen  by  the  user  consists  of  up  to  1821 
scan  lines,  each  with  3822  infrared  elements  and  15288  visible 
elements  across  the  line.   Each  infrared  element  is  referenced 
to  a  coincident  array  of  8  by  4  visible  picture  elements 
(pixels) .   The  referencing  along  each  scan  line  is  accomplished 
by  the  SDB.   The  referencing  across  scan  lines  is  governed  by 
the  VISSR  step  scan  characteristics. 

For  each  scan  line  the  actual  VISSR  scan  line  number  and  the 
value  of  "Beta"  are  contained  in  a  "documentation  block"  that 
is  prefixed  to  each  line  of  stretched  VISSR  data.   The  documen- 
tation block  also  contains  values  of  the  parameters  that  speci- 
fy the  orbit,  the  attitude,  and  the  misalignment  angles.   This 
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block  is  described  in  the  appendix  of  reference  (6).   These 
parameters  are  computed  at  the  central  facility  (VIRGS)  and 
transferred  to  the  SDB.   They  are  intended  to  furnish  the  users 
of  VISSR  data  with  the  information  needed  to  perform  predictive 
earth-location.   Clearly  any  significant  errors  in  these  para- 
meters force  the  user  to  either  work  with  mislocated 'image 
data,  or  alternatively,  to  do  the  navigation  independently.   If 
the  orbit-attitude  documentation  block  contains  precise,  or  the 
best  achievable,  values  on  a  highly  reliable  basis,  the  user 
needs  only  a  standardized  geometry  routine  to  which  the  docu- 
mented parameters  are  passed  as  inputs. 

Of  interest  to  some  types  of  users  is  the  NESS  standard  full 
disc  grids.   They  are  embedded  by  the  SDB  as  the  9th  bit  of 
each  8  bit  infrared  (IR)  sample.   The  locations  of  the  NESS 
standard  grid  points  are  generated  by  VIRGS  and  transferred  to 
the  SDB.   For  accuracies  within  the  resolution  of  an  IR  pixel 
these  grids  can  be  used.   Thus  a  user  can  avoid  the  large  com- 
putational task  of  applying  the  geometry  transform  to  some 
30,000  grid  points  that  outline  political  and  geographic  bound- 
aries.  For  users  whose  applications  require  full  visible  res- 
olution grid  point  location,  the  documented  orbit  and  attitude 
parameters  should  be  used. 

Figure  1  illustrates  the  overall  data  flow. 

D.   The  Earth-Location  Capability 

The  foundation  of  a  VISSR  earth-location  capability  is  the 
geometrical  transformation  and  its  inverse.   Given  the  earth- 
location  parameters  (the  satellite's  position  and  attitude  and 
the  VISSR1 s  misalignment  angles)  the  transformation  maps  the 
J    sample  on  the  I    scan  line  to  the  latitude  and  longi- 
tude of  the  point  on  the  earth  to  wnich  that  sample  corresponds 

The  most  convenient  and  accurate  method  for  determining  the 
values  of  the  earth-location  parameters  uses  earth  landmarks 
and  stars  that  are  discernable  in  the  VISSR  image.   The  image 
coordinates  (I,  J)  of  known  stars,  over  an  interval  of  time, 
specify  the  values  of  the  satellite's  attitude  and  the  VISSR' s 
misalignment  angles.   These  values  and  the  coordinates  (I,  J) 
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Figure  1.   Overall  data  flow  for  earth  locating  VISSR  dat 


of  one  or  more  recognizable  landmarks  (measured  during  a  time 
interval  close  to  the  the  time  interval  of  the  star  measure- 
ments) provide  a  solution  for  the  satellite's  position  as  a 
function  of  time,  i.e.,  for  its  orbit. 

Using  star  measurements  in  the  navigational  process  separates 
the  satellite's  attitude  and  misalignment  determination  from 
the  orbit  determination.   This  disjointness  greatly  simplifies, 
and  enhances  the  efficiency  of,  the  orbit/attitude 
determination  process.   These  advantages  are  maintained  even 
when  more  sophisticated  mathematical  models  are  employed.   This 
method  also  simplifies  the  structure  of  the  software  and 
streamlines  the  process  so  that  the  cpu  time  required  for  the 
orbit/attitude   computation  is  negligible.   Finally,  the  use 
of  stars  eliminates  the  need  for  the  expensive  operation  of 
ground   stations  that  provide  ranging  data  which  are  currently 
being  used  to  determine  the  orientation  of  the  orbital  plane. 

Implicit  in  the  technique  for  determining  attitude  from  star 
positions  and  orbit  from  landmark  positions,  are  mathematical 
models  of  the  dynamic  behavior  (or  motion)  of  the  spin  axis  and 
position.   Such  models  propagate  the  motion  from  initial 
(epoch)  values  for  the  parameters.   For  GOES  both  the  orbit  and 
attitude  models  are  relatively  simple  owing  to  the  well-behaved 
motion  of  the  satellite  and  its  spin  axis.   Having  relatively 
simple  models  suggests  that  if  the  orbit/attitude  determination 
methods  were  cleverly  designed,  they  could  be  implemented  as  a 
small,  efficient  software  package. 

The  less  sophisticated  models,  i.e.,  those  that  account  for 
fewer  effects,  can  only  be  used  to  accurately  predict  satellite 
motion  over  short  periods.   Real-time  data  users  need  accurate 
parameters  only  a  few  minutes  in  advance,  say  in  the  prior 
image  frame.   The  operators  of  the  central  facility,  to  avoid  a 
continual  real-time  activity  of  supporting  the  SDB ,  would  pre- 
fer a  once-per-day  cycle.   This  requires  at  least  a  24  hour 
predictive  capability.   Being  able  to  predict  over  several  days 
or  a  week  has  some  potential  advantages  in  procedural  efficien- 
cy.  It  would  support  such  things  as  maneuver  planning  and 
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batch  processing  of  the  grid  point  locations  during  off  hours. 
In  brief,  the  prediction  capability  should  be  a  minimum  of  one 
day,  and  preferably  several  days. 

Since  maneuvers  of  the  satellite  occur  frequently,  and  inevi- 
tably cause  some  degree  of  discontinuity  in  the  earth-location 
process,  an  extremely  important  aspect  of  the  earth-location 
capability  is  quick  recovery  of  full  accuracy  following  a  ma- 
neuver.  Predictions  of  the  post-maneuver  orbit  and  attitude 
are  done  but  cannot  be  expected  to  be  highly  accurate.   The 
modeling  of  maneuvers  is  not  precise  and  many  maneuvers  are  not 
conducted  as  planned.   Using  the  "stars  and  landmarks"  naviga- 
tion technique  immediately  following  a  maneuver  to  make  succes- 
sively improved  estimates  of  the  orbit  and  attitude  allows  re- 
covery of  full  accuracy  within  a  few  hours.   The  amount  of  time 
required  depends  on  the  time-of-day  that  the  maneuvers  were 
performed,  relative  to  the  time  when  the  stars  and  landmarks 
are  visible.   Orbit/attitude  determination  methods  should  en- 
able the  best  possible  estimates  of  parameters  with  a  minimal 
set  of  observations. 

Finally,  the  central  facility  must  generate  and  transfer  the 
necessary  parameters  to  the  SDB.   These  parameters,  when  their 
values  are  precise,  allow  the  SDB  to  center  the  scan  lines;  and 
by  using  these  parameters,  users  can  perform  accurate  geometric 
transforms  with  standard  geometry  routines. 
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III.   MATHEMATICAL  ALGORITHMS  AND  MODELS  FOR  THE 

EARTH  LOCATION  PROCESS 

A.  Preview 

The  algorithms  for  navigating  geosynchronous  satellite  images 
will  now  be  described.   These  algorithms  use  the  positions  of 
recognizable  stars  and  earth-based  landmarks  measured  in  VISSR 
image  frames  to  determine  the  attitude,  the  misalignment  angles 
and  the  orbit  parameters.   The  set  of  algorithms  is  divided 
into  three  parts:   (1)  the  transformation  of  image  coordinates 
to  earth  coordinates,  and  the  corresponding  inverse  transfor- 
mation, using  values  for  the  attitude,  misalignment,  and  orbit 
parameters,  (2)  the  determination  of  attitude  and  misalignment 
parameters  from  measured  correspondences  for  stars  between 
inertial  pointing  vectors  and  image  frame  positions,  and  (3) 
the  determination  of  orbit  parameters  from  the  measured  image 
positions  of  identifiable  earth-based  landmarks  and  attitude 
and  misalignment  parameters. 
B.  Geometric  Transformations  Between  Image  and  Earth  Coordinates 

The  approach  for  transforming  image  coordinates  to  earth 
coordinates  is  now  outlined.   The  details  of  the  mathematical 
formulation  of  the  transformation  follow  shortly. 

The  first  step  is  to  use  the  attitude  and  misalignment 
parameters  to  find  the  inertial  pointing  direction  of  the  spin 
scan  camera  (e.g.,  the  VISSR).   The  inertial  pointing  must  be 
found  as  a  function  of  line  number,  element  number  and  image 
start  time.   Next,  the  position  of  the  satellite  as  a  function 
of  time  is  found.   For  given  line  and  element  numbers  the 
satellites  orbital  position  and  the  spin  scan  camera's  inertial 
pointing  direction  are  known,  and  it  is  a  straight  forward 
procedure  to  determine  the  point  on  the  earth's  surface  being 
observed  by  the  satellite's  camera.   This  is  done  by  projecting 
a  ray  in  the  inertial  pointing  direction  from  the  satellite 
position  and  determining  the  point  where  it  would  intersect  an 
oblate  spheriod  (fig.  2) . 
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Figure  2. 

The  inverse  problem,  i.e.,  determining  the  line  and  element 
number  of  the  sample  which  observes  a  particular  point  on  the 
earth's  surface  in  a  given  image,  is  more  complicated.   The 
complications   arise  from  not  knowing  the  time  that  the  point 
in  question  is  being  observed.   Knowing  the  time  is  essential 
since  the  inertial  pointing  direction  of  the  vector  from  the 
satellite  to  a  particular  earth  location  varies  as  a  function 
of  time  due  to  orbital  motion  of  the  satellite  and  the  rotation 
of  the  earth  with  respect  to  inertial  space.   Since  this  time 
is  not  known,  a  guess  must  be  made.   Using  this  "first  guess" 
one  normalizes  the  inertial  vector  from  the  satellite  to  the 
point  on  the  earth's  surface  into  an  inertial  pointing  vector 
and  then  computes  the  line  and  element  number  at  which  the 
pointing  direction  of  the  spin  scan  camera  was  coincident  with 
this  pointing  vector.   In  general,  the  time  at  which  this  par- 
ticular sample  was  taken  is  different  from  the  first  guess. 
The  new  time  is  taken  and  this  procedure  is  repeated  iterative- 
ly  until  two  consecutive  times  agree  to  within  the  time  between 
two  successive  scan  line  sweeps. 

The  discussion  shows  that  three  basic  computations  must  be 
made  in  order  to  transform  image  coordinates  to  earth  coordi- 
nates and  vice  versa.   These  are:   (1)  the  determination  of  the 
inertial  pointing  direction  of  the  spin  scan  camera  as  a  func- 
tion of  line  number,  element  number,  picture  time,  beta  count, 
and  misalignment  and  attitude  parameters,  (2)  the  determination 
of  the  intersection  between  a  ray  and  the  surface  of  an  oblate 
sphere,  and  (3)  the  determination  of  the  line  and  element  num- 
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ber  at  which  the  axis  of  view  of  the  spin  scan  camera  is  paral- 
lel to  a  specific  inertial  pointing  vector  for  a  given  set  of 
beta  counts  and  misalignment  and  attitude  parameters. 

The  reader  should  be  aware  that  "spin  scan  camera"  is 
substituted  for  "VISSR"  to  connote  the  generality  of  the 
formulation..   Also  the  "pointing  direction"  of  a  spin  scan 
camera  is  the  axis  of  the  optical  field  of  view  of  the  camera 
or,  for  brevity,  the  "axis  of  view." 

Figure  3  illustrates  the  elements  of  the  basic  geometry. 
1.   Determination  of  Pointing  Direction  of  the  Spin  Scan  Camera 

The  pointing  direction  of  the  spin  scan  camera  is  found  as  a 
function  of  time,  beta  counts,  and  line  and  element  numbers. 
To  determine  this  one  must  have  an  accurate  attitude  state, 
ice.,  the  attitude  itself  and  the  precession  of  the  attitude. 
One  must  also  know  the  parameters  that  determine  the  state  of 
misalignment  between  the  satellite's  body  axis  and  its  actual 
spin  axis.   Those  parameters  establish  the  stepping  path  of  the 
spin  scan  camera  as  a  function  of  line  number. 

Nominally  the  spin  axis  and  the  body  axis  of  the  satellite 
are  expected  to  coincide.   However,  in  general,  some  small 
misalignments  are  observed.   The  misalignment  causing  a  bias  in 
the  line  positions  where  identifiable  stars  and  earth-based 
landmarks  appear  in  the  image  is  designated  as  pitch  and  param- 
eterized with  the  angle  C  .   The  misalignment  that  causes  a 
bias  in  the  element  direction  is  designated  as  roll  and  parame- 
terized with  the  angle  p  and,  finally,  the  misalignment  causing 
skew  in  the  element  direction  as  a  function  of  line  number  is 
designated  as  yaw  and  parameterized  with  the  angle  r|  . 

The  pointing  position  at  the  start  of  each  scan  line  is  ref- 
erenced to  the  sun  pulse  detection  for  that  scan  line  by  the 
VISSR  to  sun  pulse  detector  sweep  angle,  T  ,  the  current  beta 
count,  p  ,  and  whatever  element  bias  and  skew  that  are  present. 
It  is  natural  to  study  the  position  of  the  spin  scan  camera  as 
a  function  of  scan  line  number  at  those  instances  when  the  sun 
pulse  is  detected.   Two  satellite  centered  coordinate  systems 
are  created  for  this  study.   The  first  coordinate  system  has 
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its  z-axis  pointing  opposite  the  spin  axis  and  its  x-axis 
pointing  at  the  projection  of  the  sun  in  the  spin  plane--the 
plane  perpendicular  to  the  spin  axis  of  the  satellite  and  which 
passes  through  the  center  of  the  satellite  (fig.  4).   The  y- 
axis  is  added  to  form  a  right-handed,  orthogonal  coordinate 
system. 
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Figure  4 
At  the  time  the  sun  pulse  is  detected,  the  projection  of  the 
axis  of  view  of  the  sun  pulse  detector  in  the  spin  plane  coin- 
cides with  the  x-axis  of  this  coordinate  system.   Initially, 
one  may  assume  that  the  body  axis  of  the  satellite  coincides 
with  the  spin  axis  of  the  satellite.   In  this  case  the  projec- 
tion of  the  axis  of  view  of  the  spin  scan  camera  into  the  spin 
plane  coincides  with  the  vector  (cos(y) ,  sin(y) ,0)    where  y   is 
the  angle  between  the  projections  of  the  axes  of  view  of  the 
sun  pulse  detector  and  the  spin  scan  camera  in  the  spin  plane. 
Under  the  assumption  that  the  body  and  spin  axes  of  the  satel- 
lite coincide,  the  axis  of  view  of  the  spin  scan  camera,  at  the 
time  of  sun  pulse  detection,  projects  onto  the  x-y  plane  (the  spin 
plane)  at  an  angle  y    from  the  x-z  plane  of  the  coordinate  sys- 
tem.  Since  the  scan  line  number  specifies  the  position  of  the 
spin  scan  camera  above  and  below  the  x-y  plane,  one  can  now 
express  the  position  of  the  axis  of  the  spin  scan  camera  at  the 
time  of  sun  pulse  detection  as  a  function  of  scan  line  number. 
First  set  c^  to  be  the  picture  center  line,  r^  to  be  the 
radians  per  line  and  l    to  be  the  line  number.   The  pointing 
direction  (axis  of  view)  of  the  spin  scan  camera  is  given  by 
the  vector: 
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cos(y) .cos( (c^ -I) .  r^) 

sin(y)  .cos(  (c;  -  i)  .  rf  )  (2.1) 

sin(  (Cjj-1)  *r£  ) 

Since  the  body  axis  does  not,  in  general,  coincide  with  the 
spin  axis  of  the  satellite?  a  basis  for  describing  the  motion 
of  the  axis  of  view  of  the  spin  scan  camera  with  respect  to  the 
above  established  coordinate  system  is  developed  next.   Specif- 
ic misalignment  parameters  can  be  related  to  pictorial  effects 
such  as  a  line  bias,  an  element  bias  and  a  skew  bias  across  the 
elements  as  a  function  of  line  number.   To  have  this  direct 
correspondence  between  observable  effects  and  angular  param- 
eters a  second  satellite  centered  coordinate  is  defined.   The 
z-axis  again  points  opposite  the  spin  axis  and  the  x-axis  coin- 
cides with  the  (cos(y) ,  sin(y) ,  0)  vector  of  our  first  coordin- 
ate system.   Hence,  if  the  body  axis  does  in  fact  coincide  with 
the  spin  axis  of  the  satellite,  the  x-axis  of  our  second  coor- 
dinate system  coincides  with  the  axis  'of  view  of  the  spin  scan 
camera  at  the  center  line  position  when  the  sun  pulse  is 
detected.   This  coordinate  system  is  arrived  at  from  the  first 
coordinate  system  by  rotating  the  x-axis  around  the  z-axis 
towards  the  y-axis  by  the  angle  7. 

Within  this  coordinate  system  the  pointing  direction  of  the 
spin  scan  camera  as  a  function  of  scan  line  number  at  the  time 
of  sun  pulse  detection  is  as  follows: 

'  cos(  (c^-  f)  •  rp  \ 

0    J  (2.2) 

sin(  (C|  -  i )  .rf  )  i 

A  perturbation  matrix  M  enables  one  to  describe  the  pointing 
direction  (axis  of  view)  of  the  spin  scan  camera  as  a  function 
of  scan  line  number  for  the  case  of  misalignment  between  the 
spin  axis  and  body  axis  of  the  satellite.   The  misalignments 
are  so  small  that  they  can  be  broken  down  into  three  essential 
degrees  of  freedom:   (1)  rotation  around  the  y-axis  of  the 
z-axis  towards  the  x-axis  being  positive  "pitch"  (it  makes  the 
earth  appear  higher  in  the  image) ,  (2)  rotation  around  the 
x-axis  of  the  z-axis  towards  the  y-axis  being  positive  "yaw" 
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(it  skews  the  picture  making  the  top  of  the  earth  appear  shift- 
ed to  the  right  with  the  bottom  of  the  earth  appearing  shifted 
to  the  left  in  the  image)  and  (3)  rotation  around  the  z-axis  of 
the  y-axis  towards  the  x-axis  being  positive  "roll"  (it  makes 
the  earth  appear  to  move  towards  the  left  in  the  image) .   The 
order  in  which  one  corrects  for  these  misalignment  angles  is 
really  arbitrary.   However,  one  ordering  must  be  selected  to  be 
definitive.   Furthermore,  because  the  misalignment  angles  are 
so  small,  the  effect  of  applying  these  several  misalignment 
effects  in  the  wrong  order  would  be  secondary  and  undetect- 
able.  The  misalignment  matrix  is  defined  as: 

(cos(p)  sin(P)  o\  A     0       o\  /cos(£)0  sin(C)  \ 
-sin(p)  cos(P)  0.0  cos(t!)  sin(T!)  J  .(   0    1    0     I  (2.3) 
0       0     1 7  V 0  -sin(ri)cos(Ti)/  \-sin(U0  cos(U/ 
where  £,  is  pitch,  r\    is  skew  and  p  is  roll  misalignments  respec- 
tively.  Matrix  M,  follows  the  conventions  used  in  Mottershead 
and  Phillips  (1) . 

The  pointing  direction  of  the  spin  scan  camera  as  a  function 
of  scan  line  number  accounting  for  misalignment  is  as  follows: 

/cos(  {ci-l).rl  )  \ 
M  .  [        0       j  (2.3a) 

\  sin( (Cj - t)  .re ) / 
The  pointing  direction  of  the  spin  scan  camera  as  a  function  of 
scan  line  number  and  accounting  for  misalignment  in  the  first 
satellite  centered  coordinate  system  is  as  follows: 
'cos(y)   -sin(y)  0\       /cos(  (c^  -  ft)  .  r^  )  \ 
sin(y)    cos(y)  0   .  M  .  0      J        (2.4) 

0      0       1/      \sin(  (c^-H  -r^)/ 
A  way  of  envisioning  these  effects  in  the  acquired  image  data 
is  now  presented.   Superimpose  a  spherical  latitude-longitude 
grid  around  the  spacecraft  with  the  equatorial  plane  coinciding 
with  the  spin  plane  of  the  spacecraft.   Portrayed  is  the  scan- 
ned portion  of  this  grid,  +10°  to  -  10°  (actually  a  little  bit 
more) ,  on  paper  with  a  grid  of  horizontal  and  vertical  lines; 
the  horizontal  lines  represent  constant  latitudes  corresponding 
to  constant  line  numbers  and  the  vertical  lines  represent  con- 
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stant  longitude  positions.   A  nominal  rectangular  shaped  object 
within  this  framework  is  represented  in  figure  5. 

TOP  OF  PICTURE 


BOTTOM  OF  PICTURE 
(top  to  bottom  skew  about  5  ) 

Figure    5 

The  skew  from  right  to  left  with  increasing  line  number  is 
caused  by  decreasing  beta  counts  which  compensate  for  an  appar- 
ent motion  of  the  earth  from  right  to  left  due  to  orbital 
motion.   With  similar  diagrams  one  can  portray  the  effects  of 
the  pitch,  roll,  and  yaw  misalignment  parameters.   This  is  done 
by  plotting  their  effects  on  the  motion  of  the  spin  scan  camera 
as  it  steps.   The  location  of  the  sun  pulse  detector  is  desig- 
nated with  an  *  and  a  trace  is  made  of  the  path  of  the  motion 
of  the  axis  of  view  of  the  spin  scan  camera  as  the  camera  steps 
from  the  top  to  the  bottom.   The  four  plots  in  figures  6  and  7 
are  simplifications  for  a  nonspinning  spacecraft  or,  if  the 
reader  prefers,  for  the  successive  cases  where  the  sun  pulse 
detector  points  in  the  same  inertial  direction.   The  path  of 
the  axis  of  view  of  the  spin  scan  camera  is  presented  as  a 
broken  line.   For  non-zero  roll  and  yaw  misalignment  angles  the 
path  of  the  nominal  motion  of  the  camera  is  presented  with  a 
solid  line  for  comparison. 

Note  that  a  positive  pitch  has  the  camera  pointing  lower  and 
hence  makes  the  earth  appear  high  in  an  image  frame.   Positive 
yaw  has  the  camera  pointing  further  to  the  left  at  the  top  of 
the  frame.   Positive  roll  moves  the  camera  over  to  the  right 
and  hence  makes  the  earth  appear  to  be  further  to  the  left. 
The  conventions  presented  here  are  not  consistent  for  all 


18 


NOMINAL  PATH 


T 

*¥* 

NOMINAL  PATH 

OF  SPIN  SCAN  CAMERA'S 

AXIS  OF  VIEW. 


NOMINAL  ANGLE  (7)  BETWEEN 

VISSR  AND  SUN  PULSE 

DETECTOR 

Figure    6 


SUN  PULSE 
DETECTOR 


POSITIVE  PITCH 


POSITIVE  YAW 


POSITIVE  ROLL 


T 


PITCH  ANGLE 


\ 
\ 
\ 

\ 

y 

i 

< 

\ 

*  I*- 

l_tl 


YAW  ANGLE  n 


ROLL 
ANGLE  p 


Figure  7 
pieces  of  software  used  to  navigate  geostationary  satellite 
data. 

The  model  is  applied  to  find  the  pointing  direction  of  the 
spin  scan  camera  in  an  inertial  coordinate  system  from  the 
inertial  coordinates  of  the  spin  axis  vector  and  the  pointing 
vector  from  the  satellite  to  the  sun,  the  line  number,  t   ,  the 
sample  number,  s,  and  the  beta  count,  p.   The  constants  r  , 
tg ,  and  s   respectively  equal  the  radians  per  sample,  the 
time  per  beta  count  and  the  spin  rate  of  the  satellite  in 
radians  per  time.   When  sample  s  is  viewed  by  the  spin  scan 
camera  at  scan  line  2    with  beta  count  (3  ,  the  pointing  direction 
of  the  spin  scan  camera  in  the  second  satellite  centered 
coordinate  system  is  given  by 

cos(  (co  -£)  .rn ) 


cos(o )  sin(^)  0 

-sin(ff)  cos(c')  0 
0      0      1 


M 


0 


(2.5) 


sin(  (eg  -It)  .  r^  ) 
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where  a  =    (|3  .  to  .  s  )  +  (s.  r  )  -  7. 

A  rotation  matrix,  N,  is  needed  which  changes  the  representa- 
tion of  a  pointing  vector  from  the  first  satellite  centered 
coordinate  system  to  a  standard  inertial  coordinate  system, 
also  satellite  centered.   This  matrix  N  enables  one  to  find  the 
pointing  direction  of  the  spin  scan  camera  in  an  inertial  coor- 
dinate system.   Let  S_  be  the  representation  of  the  spin  axis 

p 

vector  in  the  inertial  coordinate  system,  and  let  P  represent 
the  pointing  vector  from  the  satellite  to  the  sun  as  projected 
into  the  spin  plane.   Then 

N  =  (P,  PxS  ,  -S  )  (2.6) 

Hence  the  inertial  pointing  direction  of  the  spin  scan  camera 
for  a  particular  sample  s  at  scan  line  I    with  beta  count  /3  can 
be  found  by  premultiplying  the  expression  (2.5)  with  the  matrix 
N. 
2.   Finding  Intersection  of  Pointing  Ray  and  Earth's  Surface 

Now  consider  the  problem  of  intersecting  a  ray  extended  from 
the  satellite  with  an  oblate  spheroid  representing  the  surface 
of  the  earth  (see  figure  2  again).   Assume  that  the  inertial 
coordinates  of  the  satellite  position  vector  are  (u,  v,  w)  and 
the  pointing  vector  coinciding  with  the  axis  of  view  of  the 
spin  scan  camera  for  that  particular  line  number  and  sample 
number  is  (d,  e,  f ) .   The  earth's  surface  is  represented  by  the 
equation 

(x2+y2)/a2+z2/b2  =  1 
where  a  is  the  equatorial  radius  of  the  earth  and  b  is  the 
polar  radius.   One  intersects  the  ray  (u+sd,  v+se,  w+sf)  with 
the  earth  surface  by  substituting  u+sd,  v+se  and  w+sf 
respectively  for  x,  y  and  z  in  the  equation  for  the  earth's 
surface  and  solving  for  s:   solve  the  equation 

(  (u+sd)  2+( v+se)  2)/a2  +  (w+sf)2/b2=l 
for  s.   This  equation  is  quadratic  in  s  and,  if  the  ray  does  in 
fact  intersect  this  oblate  surface,  then  the  radical  of  the 
quadratic  is  nonnegative.   In  the  case  of  the  radical  being 
positive,  the  ray  intersects  the  oblate  surface  at  two  points, 
the  one  closer  to  the  satellite  is  visible  to  the  satellite  and 
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the  point  away  from  the  satellite  is  on  the  opposite  side  of 
the  earth. 

This  completes  the  discussion  of  transforming  image 
coordinates  to  earth  coordinates.   For  the  inverse  transform, 
an  additional  technique  is  needed. 
3.   Determining  the  Line  and  Element  of  an  Inertial  Vector 

The  problem  of  converting  a  satellite  centered  inertial 
coordinate  vector  into  line  and  element  counts  for  a  particular 
satellite  image  is  now  considered.   Such  vectors  occur 
naturally  when  one  is  converting  earth  coordinates  to  image 
coordinates.   The  vector  in  this  case  is  the  normalization  of 
the  vector  extending  from  the  satellite's  position  to  the 
earth's  surface  (see  figure  2), 

The  line  and  element  numbers  at  which  the  axis  of  the  spin 
scan  camera  is  parallel  to  this  vector  is  determined  in  a  two 
step  procedure.   First  the  components  of  a  vector,  which  is 
parallel  to  this  vector,  and  has  its  base  at  the  origin  of  the 
first  satellite  centered  coordinate  system,  are  determined. 
Next  the  line  and  element  numbers  at  which  the  axis  of  the  spin 
scan  camera  coincides  with  this  second  vector  are  determined 
using  (2.5)  and  appropriate  arcsine  and  arctangent  formulae. 
The  details  for  this  transformation  are  given  in  Appendix  A. 
C.   Finding  the  Misalignment  Parameters  and  Attitude 

To  transform  image  coordinates  to  earth  coordinates,  and  vice 
versa,  the  attitude  state  and  the  misalignment  parameters  of 
the  spacecraft  are  needed.   The  attitude  state  can  be  computed 
in  three  ways?   (1)  attitude  determination  alone,  (2)  attitude 
determination  along  with  determination  of  the  pitch  misalign- 
ment parameter  and  (3)  attitude  state  along  with  a  set  of 
precession  parameters.   All  these  ways  are  useful  in  an  opera- 
tional environment.   One  may  wish  to  compute  only  the  attitude 
if  the  data  base  is  not  large  enough  to  compute  more  parame- 
ters.  For  the  same  reason  one  may  be  interested  in  computing 
only  the  attitude  and  pitch  misalignment  parameter.   Computing 
both  the  attitude  and  precession  requires  that  star  measure- 
ments be  made  over  several  days,  and  this  cannot  be  done 
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immediately  following  a  maneuver. 

The  methods  used  to  compute  the  attitude,  the  pitch  misalign- 
ment parameter,  and  the  attitude  precession  could  also  be  used 
to  compute  the  roll  and  yaw  misalignments.   However,  fitting  a 
line  through  measured  element  discrepancies,  as  a  function  of 
line  number,  and  using  the  slope  of  the  derived  line  and  the 
appropriate  intercept  as  yaw  and  roll  misalignment  angles  is  a 
simpler,  equally  effective  approach.*   This  accounts  for  all 
first  order  effects.   The  second  order  effects  are  far  below 
the  resolution  of  the  measurements.   The  effect  of  the  spin 
axis'  nutation  can  be  removed  from  the  star  measurements  by 
monitoring  the  sun  pulse  documentation  to  find  the  phase  and 
the  amplitude  of  the  sine  wave  effects  of  this  nutation.   The 
main  impact  of  this  effect  is  to  compress  and  stretch  the 
images  in  the  vertical  direction.   This  effect  should  be 
removed. 

The  data  base  for  determining  the  attitude  state  and  pitch 
misalignment  consists  of  inertial  pointing  vectors  associated 
with  image  line  and  element  pairs  at  which  there  is  an  identi- 
fiable star  or  a  recognizable  earth-based  landmark.   In  the 
case  of  an  earth-based  landmark  one  assumes  that  a  correct  set 
of  orbital  parameters  are  available  to  use  in  computing  the 
inertial  vector  from  the  satellite  to  the  earth-based  land- 
mark.  In  the  case  of  a  star   measurement,  the  inertial  point- 
ing vector  is  computed  from  the  right  ascension  and  declination 
of  the  star.   Here  for  brevity  only  star  measurements  are  con- 
sidered for  attitude  determination.   There  is  no  intention, 
however,  to  obviate  the  use  of  earth-based  landmarks  for  atti- 
tude determination. 


*Note:   Currently  on  the  VIRGS  system,  roll  and  yaw  are  com- 
puted manually.   First  the  element  discrepencies  of  two  star 
measurements  at  different  scan  lines  are  extrapolated  linearly 
to  the  element  discrepancy  at  the  center  scan  line  of  the  pic- 
ture.  This  element  discrepancy  is  then  converted  into  an  angle 
and  entered  as  roll.   Next  the  ratio  of  the  difference  of  ele- 
ment discrepencies  to  the  difference  in  scan  line  is  computed 
in  terms  of  angular  measures  and  entered  as  yaw. 
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Consider  now  the  case  of  having  several  star  measurements 
with  the  spin  axis  vector  essentially  unchanged  during  the  time 
period  in  which  the  measurements  were  taken.   Represent  the 
pointing  direction  of  the  spin  axis  with  the  unit  vector  (a,  b, 
c)  and  assume  k  star  measurements  consisting  of  the  quadruples 
(b^,    a^,     2^,    si)  with  i=l,...,k  where  6^,    a ^,     2^, 
and  s.  are  respectively  the  declination,  right  ascension,, 
line  number,  and  sample  number  of  the  star  measurement.   Then, 
by  setting  c«  and  r „    equal  to  the  picture  center  line  and 
radians  per  line,  the  following  relationship  is  obtained: 
a.cos^j)  •  cos{a  .)  +b.  cos(  L,  .)  .  sin(  q. )  +c.  sin(  £ . ) 

=cos(7r/2  +  r£  •  (Cjg  -£i)  ) 
for  i=l,...,k.   It  may  appear  that  at  least  three  measurements 
are  required  to  solve  for  the  three  entries  of  the  spin  axis 

vector,  but  in  fact  one  needs  only  two  measurements  along  with 

2   2   2 
the  constraint  that  a  +b  +c  =1.   Let  (x.,  y.,  z.) 

i  2  l    l 

for  i=l,...,k  be  the  unit  vectors  pointing  at  the  stars  in 
inertial  coordinates.   See  figure  8. 
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Figure  8 

Mathematically  the  task  is  to  minimize  the  following  expression: 

k 

S(a,  b,  c)  =  2  (ax.  +by.  +cz  -cos((^  .)  )  (3.1) 

i=l       x    1        x 

2   2   2 
subject  to  the  constraint  a  +b  +c  =1  where  tf .  equals 
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V2+r^  (c£~£-j)  an<^  where  (a,  b,  c)  is  nearly  parallel  to  the 
vector  (0,  0,  -1)  which  points  opposite  to  the  earth's  spin 
axis  to  satisfy  the  operational  constraint  of  centering  the 
earth  in  each  image.  The  details  for  minimizing  the  sum  in 
(3.1)  are  given  in  Appendix  B. 

A  nonzero  line  misalignment  parameter  C  changes  the  angle 
between  the  spin  axis  and  the  axis  of  the  spin  scan  camera  by 
the  amount  C  .   in  this  case  it  is  appropriate  to  find  (a,  b,  c) 
and  C  to  minimize: 

k  2 

S(a,  b,  c,  C)=  2  (axi+byi+czi-cos(^i~C)  )        (3<,2) 

i=l 

2   2   2 
sub3ect  to  the  constraint  a  +b  +c  =1.   The  details  of  the 

mathematics  leading  to  such  a  minimizatin  are  given  in  Appendix 

B. 

The  element  and  skew  misalignment  parameters  can  be  ignored 
because  they  have  a  negligible,  secondary  impact  on  the  angle 
between  the  spin  axis  vector  and  the  spin  scan  camera. 

The  next  problem  is  to  find  the  precession  of  the  satellite. 
Since  this  is  a  study  of  parameters  varying  with  time,  it  is 
assumed  that  the  star  measurements  have  additionally  a  record 
of  time,  t.  associated  with  the  i    position.   An  under- 
lying assumption  here  is  that  the  precession  can  be  modeled 
well  by  constraining  the  projection  of  the  spin  axis  vector 
into  the  equatorial  plane  to  move  in  a  straight  line  with  the 
rate  of  this  motion  being  linear  in  time.   This  accounts  for 
the  first  order  effects  of  the  precession.   Figure  9  shows  the 
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model:   Here  the  (p,r)  position  is  the  projection  of  the  spin 

axis  vector  at  time  equal  to  zero,  and  the  (q,s)  vector  gives 

the  direction  of  the  precession  movement  with  the  length  of  the 

vector  corresponding  to  the  precession  rate. 

The  normalization  constraint  is  that  the  spin  axis  vector 

point  at  (p+qt,  r+st,  -  .(1-  (  (p+qt)  +  (r  +  st)  )  )A)    for  all  times  t. 

Consider  then  the  problem  of  minimizing  the  following  function: 

S(p,q,r,s)=  (3.3) 

k  1/ 

S  (  (p+qt.)  x.  +  (r+st.)  y.-(l-(  (P+qt.)  2+(r+st.)  2)  )  2  z.-cos^.)  )  ? 

The  mathematics  of  this  minimization  are  similar  to  that  pre- 
sented in  Appendix  B.   If  one  so  desires,  the  pitch  misalign- 
ment can  be  included  in  these  calculations  with  the  main  effect 
being  that  a  five  by  five  matrix  instead  of  a  four  by  four 
matrix  must  be  inverted  in  the  minimization  process. 

At  least  two  days  of  measurements  are  required  to  apply  this 
scheme  and  three  days  of  data  would  certainly  provide  more 
confidence  in  one's  results. 

D.   Orbit  Determination  Algorithms 

Next  to  be  presented  are  algorithms  to  determine  the  orbit  of 
a  geosynchronous  satellite  by  using  position  measurements  in 
image  condinates  of  earth-based  landmarks  and  assuming  one  al- 
ready has  the  Satellite's  attitude  and  misalignment  parameters 
(i.e.,  from  star  positions).   To  get  started,  the  geometry  and 
algorithms  will  be  illustrated  briefly  with  diagrams  and  dis- 
cussions that  emphasize  the  simplicity  of  the  approach  while 
leaving  the  technical  details  for  later  formulation. 

By  knowing  the  attitude  and  misalignment  angles  of  the  space- 
craft one  can  determine  the  inertial  pointing  direction  for  any 
line-element  sample.   Hence,  when  a  landmark  is  observed  on  the 
earth's  surface,  the  direction  from  the  satellite  to  the  land- 
mark is  known.   Because  the  direction  from  the  landmark  to  the 
satellite  is  opposite,  one  can  also  determine  the  direction 
from  the  landmark  to  the  satellite  by  changing  the  sign  of  the 
components  of  the  pointing  vector.   Hence  a  ray  locus  on  which 
the  satellite  must  lie  extends  from  the  landmark  towards  the 

25 


satellite.   By  intersecting  this  ray  (see  fig.  10)  with  a 
sphere  whose  radius  nominally  equals  the  distance  of  a  geosyn- 
chronous satellite  from  the  earth's  center  and  with  its  center 
coinciding  with  the  earth's  center,  a  vector  approximating  the 
satellite's  position  is  found. 
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Figure    10 

Thus  far  the  height  of  the  satellite  has  not  been  deter- 
mined.  Only  a  nominal  height  has  been  assigned  so  as  to  enable 
the  determination  the  subsatellite  point.   If  the  landmark 
measurement  being  used  is  close  to  the  subsatellite  point,  the 
method  of  determining  the  subpoint  position  by  intersecting  a 
sphere  with  a  ray  introduces  almost  no  error.   The  goal  is  to 
track  the  motion  of  the  subpoint  in  order  to  determine  a  set  of 
Keplerian  orbit  parameters. 

The  task  of  finding  Keplerian  orbit  parameters  is  split  into 
two  parts:   (1)  determining  the  plane  in  inertial  space  con- 
taining the  orbital  track  and  (2)  approximating  the  angular 
motion  of  the  satellite  around  the  center  of  the  earth  as  a 
function  of  time.   For  each  landmark  measurement,  one  is  able 
to  approximate  a  pointing  vector  in  inertial  coordinates  which 
points  from  the  center  of  the  earth  towards  the  satellite  posi- 
tion.  Index  these  pointing  vectors  and  their  associated  times 
by  i  yeilding  (t.,  x.,  y.,  z.)  for  i  equal  to  from  one 
to  the  number  of  landmark  measurements  being  considered. 

To  determine  the  orbit  plane  one  determines  a  vector  that  is 
perpendicular  or  nearly  perpendicular  to  all  these  inertial 
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vectors.   Since  the  satellite  moves  along  with  the  earth's  mo- 
tion from  west  to  east,  the  right-hand  rule  determines  that  the 
orbit  plane  vector  points  nearly  parallel  to  the  earth's  spin 
axis.   As  with  the  attitude  determination,  let  (a,  b,  c)  be  the 
inertial  unit  vector  that  is  perpendicular  to  the  orbital 
plane.   Hence  this  vector  must  satisfy  the  constraints 

ax.  +  by.  +  cz.  =  0  for  i=l,  . . . ,N 
1     2  i      i 

where  N  is  the  number  of  landmark  measurements  under  considera- 
tion.  Because  of  measurement  errors  and  other  imprecisions  in 
the  system  it  is  appropriate  to  minimize  the  following  sum: 


N 

S(a,b,c)  =z  (ax.  +  by.  +  cz.  )' 
.  ,11      i 
i=l 


(4.1) 


The  mathematical  form  of  this  problem  is  identical  to  that 
presented  with  equation  (3.1) .   The  methods  used  to  minimize 
this  expression  are  similar  to  methods  used  to  minimize  (3.1). 
Referred  to  Appendix  B  for  further  information. 

Figure  11  illustrates  the  orbit  plane  perpendicular  that  is 
determined  by  these  calculations. 
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Figure  11 
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The  orbit  plane  perpendicular  (a,b,c)  determines  the  plane  in 
which  the  satellite  moves.   However,  to  obtain  the  more  famil- 
iar Keplerian  orbit  parameters,  one  must  manipulate  the  compo- 
nents of  this  unit  vector  to  obtain  the  ascending  node  and 
inclination  of  the  orbit.   By  definition  the  pointing  vector 
(a,b,c)  is  parallel  to  the  vector 

(  sin(£2)  ■>  sin(  i)  ,  -sin(  i)  «  cos(ft)  ,  cos(i)  )  (4.2) 

where  fi  is  the  ascending  node  and  i  is  the  inclination  of  the 
orbit.   Since  the  inclination  is  always  taken  to  be  positive, 
the  ascending  node  and  the  inclination  can  be  determined  as 
follows  (in  Fortran)  : 

Q  =  ATAN2(a,  -b)  and  (4.3) 

i  =  ASIN(SQRT(a2+b2)  )  . 
This  completes  the  discussion  of  determining  the  orbit  plane 
for  the  satellite. 

In  order  to  study  the  angular  motion  of  the  satellite  within 
its  orbit  plane,  a  planar  coordinate  system  is  set  up  in  the 
orbital  plane.   The  x-axis  points  from  the  center  of  the  earth 
to  the  ascending  node  and  the  y-axis,  also  lying  in  the  orbital 
plane,  is  90°  beyond  the  x-axis  in  the  direction  of  satellite 
motion.   Refer  to  figure  12. 

Angles  in  this  planar  coordinate  system  are  measured  as  going 
from  the  x-axis  towards  the  y-axis,  i.e.,  in  the  direction  of 
satellite  motion.   Specify  angular  positions  occuring  at  time 
t.  as  a  .  .   The  argument  of  perigee  is  specified  by  oj,time 
is  specified  by  t,  and  the  time  when  the  satellite  is  at 
perigee  is  specified  by  t  .   Within  this  framework  the  motion 

XT 

of  a  satellite  in  a  Keplerian  orbit  is  well  represented  by  the 

following  equation: 

a  -oj-2  e .  sin(ff)  .  cos(uj)  +  2e  .  cos  {a)  .  sin  (w)  =>\-it-t    ),     (4.4) 

where  n  is  the  mean  motion  constant.   This  representation  is 

-3 
accurate  to  within  .12  km  for  eccentricities  less  than  10 

Since  only  the  form  of  this  equation  is  of  primary  interest, 

the  derivation  verifying  the  accuracy  is  presented  in  Appendix 

C.   Normally  eq.  (4.4)  is  written  in  the  following  equivalent 

form 
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EARTH 


Figure  12 

t.  -  c,  +  c~*a.    +  c,sin(  a.)  +  c.cos(a.)         (4.5) 
1     1     2   l     3      l      4      i 

where  the  index  i  indicates  the  times  and  angular  positions 
corresponding  to  the  landmarks  and  c,=tp-u/n,  c~  =  \/x\, 
c3=-2e*cos(co)  /n,  and  c.  =  2  £•  sin(oj)  /n.   It  is  vital  that  the 
constants  |c|  be  related  directly  to  Keplerian  orbit  param- 
eters so  that  the  reason  for  determining  these  constants 
becomes  clear. 

Observe  the  following  relations: 


n  =  1/c 


2  ' 


(4.6) 


where  n  is  the  mean  motion  constant; 

.2/3/  2/3  ..    -, 

a  =  r*  k    /n  '  ,  (4.7) 

where  a  is  the  semimajor  axis,  r  is  the  radius  of  the 

earth,  and  k  is  the  gravitational  constant  of  the  earth; 

2      2  Vs 
e  =  (c3   +  c4  )/V(2c2)  and  (4.8) 

w=  ATAN2(c4,-c3) , 

where  e  and  u>  are  the  eccentricity  and  argument  of  perigee, 

respectively.   The  inverse  relationship  for  c2,c,  and  cA 

are  given  later  in  (4.13). 

The  inverse  relationship  between  the  mean  motion  constant  n 
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and  c~  is  due  to  n  being  the  multiplicative  constant  relating 
time  to  angular  position.   The  relationship  for  determining  the 
semimajor  axis  comes  from  Escabol's  Methods  of  Orbit  Determina- 
tion, henceforth  referred  to  as  EscaboL 

The  argument  of  perigee  is  determined  in  the  above  manner 
because  (c.,    -c-J  points  parallel  to  the  vector  (sin(w) , 
cos(w) )  due  to  the  eccentricity  and  mean  motion  constant  both 
being  positive.   The  formula  for  eccentricity  comes  from  the 
observation  that 

(c32+c42)/c2  =(-2€Cos(co)  )  2  +  (2esin(w))2  =  4e2.     (4.9) 
The  computation  for  finding  a  mean  anamoly  for  a  given  epoch 
is  more  involved.   First  one  finds  the  time  of  perifocal 
passage  by  substituting  w  in  (4.5)  which  is  c,  +  c?  oo + 
c^sin(w)  +  c4cos(oj)=t  .   Since  the  vector  (c3,c4)  is 
parallel  to  (-cos(w)  ,sin(w))  and  hence  perpendicular  to  (sin(oo) 
,  cos(w)) ,  the  last  two  terms  cancel,  giving  t  ^t+Cow. 

P     J.    £• 

Now,  by  definition,  the  mean  anamoly  at  an  epoch  time  t 

equals  n»(t  -t  ) .   Hence,  adopting  the  programming  symbol 
e   p 

MEANOM  for  the  mean  anamoly, 

MEANOM  =  n(t   -  C,  -  c2«w)  =  (t   -  C^) /c2  -co.      (4.10) 
This  completes  the  process  of  transforming  the  (cf  into  orbital 
elements. 

Now  consider  the  problem  of  finding  the  jcj-  from  a  set  of 
measurements  (a-,  t .) .   For  this  purpose  linear  regression 
is  applied.   The  following  expression  is  minimized: 

S=  2  (ci+c2cyi  +  c3sin(ai)  +  c4cos(oi)-t.)  .      (4.11) 
i=l 
To  solve  for  the  |c[,  set  the  following  four  partial 

derivatives  to  zero: 

"8c,(cl'  C2'    c3f    C4}  =  °' 

9c2(cl'  c2'    c3'    c4}  =  °' 

—  (clf  c2,  c3,  c4)  =  0   and 
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9S  .  \  _  n 

8c"4(Cl'  c2'  C3'  C4)    U' 

This  results  in  a  matrix  equation  for  the  c.'s  which  is  now 

presented  with  the  following  abuse  of  notation:   SM  stands  for 

sum  over  i  for   1=1,.=,  N,  A  for  the  angles  a.,    SN  for 

sin(tt.) ,  CS  for  cos(ff-) ,  T  for  the  time  t.,  and  2  means 
ii  i 

that  the  immediately  proceeding  term  is  squared. 

N       SMA     SMSN    SMCS  \   /c-\\      /SMT   \ 
SMA     SMA2    SMASN   SMACS  \  /  C2  \_ /  SMAT  \ 
SMSN    SMASN   SMSN2   SMSNCS  I*  I  c3  /  I  SMSNT  I 
SMCS    SMACS   SMSNCS  SMCS2  /  \c^/        \SMCST/ 
The  details  of  the  specific  mathematical  technique  used  to  in- 
vert this  matrix  are  given  in  Appendix  D.   Notice  here  that  the 
inverse  of  the  matrix  is  ill-conditioned  over  short  time  inter- 
vals.  This  results  from  the  inability  to  resolve  the  differ- 
ence between  a  straight  line  and  part  of  a  sine  wave  over  a 
short  arc.   For  this  reason,  to  achieve  results  comparable  in 
accuracy  to  one's  measurements,  one  should  make  landmark  meas- 
urements distributed  over  a  time  interval  at  least  18  hours 
long.   A  convenient  way  to  circumvent  making  such  a  long 
stretch  of  landmark  measurements  is  to  make  8-hour  stretches  of 
landmark  measurements  24  hours  apart. 

In  daily  operations  it  is  often  the  case  that  one  uses  the 
semimajor  axis  computed  from  yesterday's  and  today's  data  but 
only  uses  today's  data  to  compute  the  three  along-track  orbit 
parameters:   mean  anamoly?  eccentricity,  and  argument  of  peri- 
gee.  This  is  a  two-step  procedure.   First  compute  the  constant 
c~  by  inverting  the  method  used  to  compute  the  semi-major 
axis  from  the  constant  c~.   Then  compute  the  remaining  c.'s 
by  applying  linear  regression  in  the  same  manner  used  to  find 
all  four  c.'s.   In  this  case  the  resultant  three-by-three 
matrix  is  inverted  by  applying  Cramer's  Rule.   The  along-track 
orbit  parameters  are  then  found  in  the  same  manner  as  before. 
This  particular  approach  is  potentially  valuable  for  the  time 
period  immediately  following  orbit  maneuvers  if  one  can  obtain 


31 


an  estimate  of  the  orbit's  semimajor  axis  from  energy  consider- 
ations. 

The  estimate  for  the  orbit's  mean  anamoly  can  be  refined  when 
there  are  good  estimates  of  the  other  orbit  parameters  by 
requesting  that  only  one  further  parameter  be  generated.   In 
this  case  c~,    c^  and  c.    are  generated  by  the  following 
formulae: 

c2  =  k.(r/a)3/2 

c^  =  -2 • e •cacos(w)  and  (4.13) 

c 4  =  2«  e  »cisin  (w)  . 
Finally,  if  only  two  along-track  orbital  elements  are  wanted, 
linear  regression  is  used  to  solve  for  the  parameters  c,  and 
c2r  and  these  parameters  are  transformed  into  a  semimajor 
axis  and  mean  anamoly  while  at  the  same  time  making  the  eccen- 
tricity and  argument  of  perigee  equal  to  zero.   There  are  also 
ways  of  computing  c,  and  c2  while  retaining  nonzero  values 
for  the  eccentricity  and  argument  of  perigee.   These  are  cur- 
rently not  implemented. 

This  ends  the  discussion  of  orbit  determination. 
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IV.   DESCRIPTION  OF  IMPLEMENTATION 
A.   The  Orbit  and  Attitude  Determination  Software 

The  functions  of  the  software  implementation  of  the  star  nav- 
igation can  be  split  into  three  categories:   (1)  the  scheduled 
ingest  of  images  containing  either  stars  or  earth-based  land- 
marks and  the  measurements  of  the  image  locations  of  these  fea- 
tures, (2)  the  maintenance  of  navigational  files  containing  the 
measurements  and  parameters  necessary  for  navigation  and  (3) 
the  execution  of  software  that  determines  the  necessary  orbit, 
attitude  and  misalignment  parameters.   Further  details  of  the 
software  execution  for  tasks  (1)  and  (2)  are  omitted  because 
they  are  not  germane  to  the  purpose  of  this  paper.   Instead  the 
mathematical  aspects  of  measurement  errors  in  task  (2)  affec- 
ting the  scheduling  in  task  (1)  are  discussed.   These  aspects 
affect  the  scheduling  of  stars  for  the  determination  of  the 
attitude  state  and  misalignment  parameters  and  the  scheduling 
of  earth  based  landmarks  for  the  determination  of  orbital  state. 

Star  measurements  should  be  distributed  with  respect  to  time 
and  the  right  ascension  of  the  stars1  coordinates  to  provide  an 
adequate  data  base  for  generating  the  attitude  parameters,  the 
precession  parameters  and  the  misalignment  parameters.   Simi- 
larly landmark  measurements  should  be  distributed  appropriately 
with  respect  to  time  to  provide  an  adequate  data  base  for  gen- 
erating orbit  parameters.   Such  distributions  of  measurements 
are  essential  to  prevent  an  ill-conditioned  inverse  of  a  matrix 
from  blowing  up  the  impact  of  measurement  errors.   It  is  impor- 
tant to  know  what  distribution  of  particular  measurements  are 
necessary  for  the  accurate  determination  of  specific  subsets  of 
these  parameters.   Such  knowledge  enables  one  to  establish  a 
schedule  for  ingesting  stars  and  earth  based  landmarks. 

For  attitude  determination  alone,  assuming  the  precession 
rate  is  small,  at  least  two  star  measurements  are  necessary 
with  the  right  ascension  separation,  mod  (180°),  being  15°  or 
more.   Mathematically  the  ideal  case  is  to  have  the  right  as- 
censions separated  by  90°.   For  attitude  and  pitch  alignment 
determination,  three  star  measurements  are  necessary,  two  with 
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right  ascension  separation,  mod  (180°) ,  of  15°  or  more  and  a 
different  pair  (i.e.,  a  member  of  the  first  pair  matched  up 
with  the  third  star  measurement)  with  right  ascension  measure- 
ment separation,  mod  (360°) ,  of  15°  or  more.   For  determining 
the  attitude  and  precession  parameters,  four  star  measurements 
consisting  of  two  disjoint  pairs  with  each  pair  of  stars  having 
right  ascension  separation,  mod(180°) ,  of  15°  or  more  and  with 
time  separation  from  the  first  pair  to  the  second  pair  of  at 
least  18  hours.   In  this  case  the  same  stars  can  occur  in  each 
pair.   For  determining  the  roll  and  yaw  misalignment  parame- 
ters, two  star  measurements  are  needed  with  line  separation  of 
5,000  to  10,000  lines. 

Questions  have  been  raised  about  the  viability  of  routinely 
distinguishing  stars  from  background  noise.   Such  doubts  need 
to  be  examined  because  of  the  central  role  star  measurements 
play  in  a  stand-alone  navigational  system.   Since  star  measure- 
ments are  not  currently  being  used  to  determine  the  satellite's 
attitude,  an  auxiliary  system  of  ground  stations  and  a  second 
computer  has  to  be  maintained. 

It  is  vitally  important  to  consider  each  star  measurement  as 
a  single  point  in  a  consistent  data  set.   Regarding  each  star 
measurement  as  a  separate  entity  makes  the  measurement  task 
more  difficult  and  less  fruitful.   Simply  stated,  at  a  single 
time  the  spin  axis  vector  points  in  only  one  direction.   Hence, 
if  a  set  of  star  measurements  yields  large  residuals  after  a 
spin  axis  determination,  at  least  one  measurement  is  wrong. 
The  measurement ( s)  whose  residuals  are  noticely  distinct  from 
the  residuals  of  other  measuremenets  should  be  deleted.   This 
idea  can  also  be  used  predictively  since  the  change  of  attitude 
over  a  one  day  period  is  small.   Hence  the  measured  position  of 
a  star  should  be  close  to  its  predicted  position.   Use  of  this 
idea  in  a  procedual  manner  should  greatly  reduce  problems  of 
distinguishing  stars  from  noise. 

In  addition,  a  star  often  occurs  more  than  once  on  a  given 
day.   By  adding  the  image  frames  containing  a  recurrent  star, 
with  an  appropriate  element  lag  displacement,  the  signal  to 
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noise  ratio  of  the  star's  occurrence  will  be  significantly 
enhanced.   Also,  the  full  precision  (6  bit)  visible  image 
should  be  used  in  scanning  for  stars. 

For  the  determination  of  the  orbit's  inclination  and  ascend- 
ing node,  at  least  two  landmarks  are  necessary  with  time  sepa- 
ration of  4  to  6  hours.   For  the  determination  of  the  orbit's 
mean  anamoly  once  the  other  five  orbit  parameters  are  known, 
only  one  landmark  measurement  is  needed.   For  the  determination 
of  the  orbit's  mean  anamoly,  eccentricity,  and  argument  of  per- 
igee once  the  other  three  orbital  parameters  are  known,  at 
least  three  landmark  measurements  distributed  over  a  six  to 
eight  hour  interval  are  necessary.   In  order  to  determine  all 
the  orbital  parameters  at  once,  at  least  four  landmark  measure- 
ments distributed  over  a  sixteen  to  eighteen  hour  interval  are 
necessary . 

The  computer  memory  space  necessary  for  both  the  attitude  and 
orbit  determination  programs  is  less  than  16,000  words  and 
their  execution  times  are  below  1  second.   Both  programs  are 
easily  transferable  to  other  computers  with  Fortran  compilers. 
To  achieve  such  a  transfer  one  has  to  provide  files  for  a  navi- 
gational data  base,  routines  to  monitor  and  maintain  this  data 
base  and  interface  programs  between  the  navigational  data  base 
and  the  orbit  and  attitude  determination  software. 

Fortran  listings  of  the  orbit  and  attitude  determination  rou- 
tines, UPGORB  and  FINDAT,  are  provided  in  Appendix  E. 
B.  Standard  Geometry  Routine  for  Users 

Direct  readout  users  have  access  to  the  orbit,  attitude,  and 
misalignment  parameters.   These  parameters  as  determined  by 
VIRGS  are  updated  in  the  stretched  VISSR  orbit/attitude  block. 
A  description  of  the  format  and  definitions  of  this  block  is 
given  in  reference  (6). 

A  geometric  transformation  that  uses  these  parameters  as  in- 
put is  available  from  NESS  in  a  standard  form  as  two  Fortran 
routines.   "IMGLOC"  transforms  from  earth  latitude  and  longi- 
tude to  image  frame  line  and  element.   "EARLOC"  is  the  inverse 
routine,  transforming  from  VISSR  line  and  element  to  earth  lat- 

35 


itude  and  longitude.   Listing  of  these  two  routines  are  pro- 
vided in  Appendix  F. 

These  routines  do  not  use  "Beta",  the  angle  subtended  at  the 
satellite  by  the  sun  and  earth,  even  though  it  is  provided  in 
the  documentation  block.   The  SDB  uses  a  "Beta"  function  that 
is  precisely  derived  from  the  documented  orbit,  attitude,  and 
misalignment  parameters.   Therefore,  centering  of  scan  lines  on 
the  earth  disc  has  to  be  assumed  by  the  user.   This  simplifica- 
tion's benefit  for  the  user  is  that  fewer  parameters  have  to  be 
stored.   Modification  of  the  routines  to  account  for  earth  disc 
offsets  is  a  small,  straightforward  task. 

An  assembly  language  version  of  "IMGLOC"  was  implemented  on 
the  control  computer  of  the  Central  Data  Distribution  Facility 
operated  by  NOAA/NESS  for  the  GOES-TAP  service.   Using  the 
earth-location  parameters  in  the  VISSR  data  stream,  the  routine 
locates  the  image  sectors  in  near  real-time.   The  routine  has 
been  operating  since  June,  1979. 

C.   The  VISSR  Image  Registration  and  Gridding  System  (VIRGS) 

VIRGS  is  the  central  earth  location  facility.   VIRGS,  in 
terms  of  hardware  and  basic  software,  is  a  copy  of  the  McIDAS 
II  (Man-computer  Interative  Data  Access  System) .   McIDAS  II  was 
designed  and  built  by  the  Space  Sciences  and  Engineering  Center 
(SSEC)  of  the  University  of  Wisconsin  as  an  upgraded  version  of 
the  original  McIDAS  (3,4).   The  VIRGS  was  installed  by  SSEC  at 
the  World  Weather  Building  in  Camp  Springs,  Md  in  July  1978. 
The  operational  use  of  VIRGS  began  in  May,  1979. 

The  orbit/attitude  determination  techniques  and  software 
package  described  in  this  paper  were  developed  by  SPAAM,  Incor- 
porated under  contract  to  NESS  during  the  past  2  years  as  a 
follow-on  to  work  began  by  Dennis  Phillips  while  at  SSEC. 

The  VIRGS  (or  Mcldas  II)  consists  of  a  small  computer  (the 
Harris/6  with  64,000  words  of  memory),  an  on-line  40  Mbyte 
disc,  and  an  interactive  display  console  for  the  analyst.   The 
important  characteristic  of  Mcldas  is  the  elaborate  set  of 
"key-ins"  the  analyst  has  available  for  processing  and  dis- 
play.  The  demonstrations  of  accurate  VISSR  image  location  were 
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first  performed  on  the  Mcldas  in  1974  using  the  interactive 
key-ins  to  iteratively  adjust  orbit-attitude  parameters  to 
achieve  accurate  locations  of  visible  landmarks  in  the  images 
(2,5) . 

The  VIRGS  performs  the  central  earth  location  process  on  a 
daily  basis  for  each  GOES  in  four  basic  steps: 

i.  measurement  of  star  and  landmark  positions  in  the  vis- 
ible imagery  (and  sometimes  in  the  infra-red  imagery) . 
ii.  determination  of  the  misalignment  angles  and  attitude 
from  the  star  observations  and  the  orbit  from  the 
landmark  observations 
iii.  production  and  transfer  (to  the  SDB)  of  the  "Beta" 

function,  earth-location  parameters  for  the  documenta- 
tion block,  and  the  locations  of  points  in  the  NESS 
standard  grid  table,  all  for  the  upcoming  24-hour 
period, 
iv.  monitoring  of  accuracy  of  the  earth-location  through- 
out the  day,  with  updating  if  necessary. 
The  primary  roles  of  the  human  operators-analysts  in  performing 
this  daily  operation  are: 

1.  Scheduling  the  ingest  of  image  sectors  containing 
stars  and  then  the  measurement,  using  a  joystick  con- 
trolled cursor,  of  the  bright  spot  corresponding  to 
each  star. 

2.  Selecting  cloud-free  areas  on  earth  that  contain  fa- 
vored landmarks,  scheduling  the  ingest  of  the  appro- 
priate sectors,  and  then  measuring  with  the  cursor  the 
image  positions  of  the  recognizeable  features. 

3.  Judging  a  sufficient  set  of  star  and  landmark  observa- 
tions and  entering,  in  proper  sequence,  the  several 
key-ins  that  determine  misalignment  parameters,  atti- 
tude, and  orbit. 

4.  Verifying  that  error  residuals  as  output  by  various 
routines  are  suitably  small  and  that  the  various  files 
have  been  updated. 

5.  Entering  key-ins  that  execute  routines  that  produce 
files  for  the  SDB,  quality  checking  these  files,  and 
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transferring  the  files  to  the  SDB. 

6.  Routinely  checking  that  the  parameters  and  embedded 
grids  that  are  currently  on  line  (i.e.,  in  the  data 
stream)  are  suitably  accurate,  generating  updates  for 
the  SDB  as  required. 

7.  Coordinating  with  the  satellite  control  center,  on  the 
scheduling  of  maneuvers  and  special  VISSR  imaging 
sequences,  and  appropriately  changing  the  standard 
operational  schedule  to  best  deal  with  such  events. 

Clearly  the  human  plays  a  critical  role  in  the  VIRGS  opera- 
tion.  While  the  computational  tasks  are  fully  automated  as 
Fortran  routines  executable  by  simple  key-ins,  the  human  is 
required  to  supervise  their  execution,  as  a  safeguard  against 
erroneous  measurements  under  normal  circumstances  and  as  a  re- 
serve of  judgement  against  abnormal  circumstances.   Given  the 
number  of  contingencies  that  exist  daily,  the  VIRGS  implementa- 
tion appears  to  have  an  effective  mix  of  human  and  computer 
talents. 

An  important  feature  of  the  VIRGS  implementation  is  the  clos- 
ed loop  formed  when  the  VISSR  data,  appended  with  the  VIRGS 
produced  earth-location  information  (parameters  and  embedded 
grids),  are  ingested  and  displayed  at  VIRGS.   Figure  13  por- 
trays the  loop.   When  ingesting  VISSR  data,  VIRGS  is  emulating 
any  other  direct  readout  user  and  hence  has  the  capability  to 
monitor  the  performance  of  the  system,  end  to  end.   The  VIRGS 
operator  should  be  the  first  data  user  to  become  aware  of  deg- 
radations in  the  earth-location  accuracy.   In  most  cases,  the 
operator  will  be  able  to  correct  them  quickly  through  coordina- 
tion with  SDB  operators  or  satellite  control  center. 

The  issue  of  optimal  versus  minimal  sets  of  star  and  landmark 
position  measurements  is  discussed  in  section  IV-A  as  a  part  of 
the  orbit/attitude  determination  software  descriptions.   The 
screening  of  observations,  based  on  the  confidence  factor  re- 
sulting from  any  sort  of  difficulty  in  making  the  measurements, 
is  an  intuitive  manual  task  for  the  analyst.   Although  most 
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situations  are  clear-cut,  some  require  the  collection  of  addi- 
tional measurements  for  a  reliable  solution. 

D.   Rationale  for  Centralized  Service 

Since  the  earth-location  parameters  generated  by  VIRGS  at  the 
central  facility  are  appended  to  the  stretched  VISSR  data, 
other  direct  readout  users  do  not  need  to  exercise  their  own 
capability  to  generate  these  parameters,  _if  the  centrally  gen- 
erated parameters  are  sufficiently  accurate.   A  centralized 
service  is  cost  effective  since  all  data  user  can  avoid  the 
expense  of  developing  and  operating  their  own  earth-location 
capability.   The  operators  of  VIRGS  strive  for  the  greatest 
possible  accuracy  and  use  the  built-in  safeguards  as  required 
to  guarantee  a  high  reliability. 

The  users,  knowing  that  the  earth-location  parameters  are 
readily  available  in  the  data  stream,  need  only  to  employ  a 
standard  geometery  transformation  routine.   The  time  consuming 
processing  to  derive  orbit,  attitude,  and  misalignment  parame- 
ters is  performed  for  users  by  VIRGS. 

Doubt  that  1  pixel  accuracy  is  achievable  in  an  operational 
environment  exists  because  the  demonstration  that  such  accuracy 
is  achievable  has  only  been  done  on  the  McIDAS  System  at  the 
University  of  Wisconsin  and  hence  has  occured  only  in  a  univer- 
sity environment.   The  skills  to  achieve  this  accuracy  are  em- 
bodied in  one  person,  Mr.  John  T.  Young,  and  these  skills  in 
their  totality  have  not  been  transferred  to  any  other  in- 
stallation.  However,  since  the  mathematical  algorithms  cur- 
rently available  on  the  VIRGS  System  for  determining  the  satel- 
lite's orbit/attitude  states  and  misalignment  parameters  are  a 
refined  version  of  the  orbit/attitude  model  of  the  McIDAS  used 
by  Mr.  Young,  one  would  expect  that,  after  software  and  proced- 
ual  deficiencies  are  recognized  and  removed,  1  pixel  accuracy 
will  be  achieved.   The  actual  results  discussed  in  Section 
IV-E,  obtained  on  VIRGS  by  Larry  Hambrick,  lend  strength  to 
this  expectation.   One  advantage  of  having  these  algorithms  in 
the  form  of  computer  programs  (as  they  are  now  available)  is 
that  an  operator  can  objectively  generate  parameters  describing 
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satellite's  orbit/attitude  state  without  having  a  deep  under- 
standing of  orbit/attitude  relationships. 

An  accuracy  of  1  pixel,  as  used  herein,  is  defined  as  the 
root-mean-square  of  the  earth  location  errors,  line  and  element 
errors  taken  separately,  distributed  over  the  earth  disc 
throughout  a  24  hour  period.   One  pixel  is  being  used  here  as 
the  angular  distance  between  successive  pixels  in  a  full  reso- 
lution picture.   With  this  definition  one  is  not  guaranteed 
that  any  given  point  on  any  given  frame  will  be  located  to 
within  1  pixel.   A  specific  point  at  a  given  time  might  be  in 
error  by  several  pixels.   However,  loosening  this  requirement 
further,  say  to  two  or  more  pixels  rms,  significantly  increases 
the  possibility  of  larger  errors.   In  terms  of  ground  distance, 
a  one  pixel  error  spans  one  kilometer  at  the  subpoint  but  spans 
two  kilometers  50°  from  the  subpoint.   Finally,  it  is  noted 
that,  although  1  pixel  accuracy  rms  is  probably  acceptable  to 
all  current  users,  even  sharper  accuracies  are  possible  with 
the  precision  alignment  capabilities  inherent  in  the  data 
stream. 

The  software  currently  available  on  the  VIRGS  System  has  been 
verified  to  approach  1  pixel  accuracy.   Even  though  routine 
operators  are  not  now  expected  to  attain  this  level  of  accura- 
cy, it  appears  likely  that  information  transfer  and  training 
will  improve  the  routine  operation.   The  goal  is  to  make  the 
routinely  achieved  accuracy  a  characteristic  of  the  software 
package  alone.   This  can  be  accomplished  by  improving  the  soft- 
ware packages  to  reduce  the  reliance  on  an  operators' s  judge- 
ment and  by  eliminating  the  impact  of  an  operator's  bias 
through  appropriate  training. 

E.   An  Error  Budget  and  Some  Actual  Results 
1.   Error  Budget 

The  accuracy  goal,  stated  here  for  overall  system  perform- 
ance, is  one  full  resolution  visible  pixel,  which  is  a  21  jara- 
dian  angular  interval  at  the  satellite  that  translates  into  0.8 
km  at  the  subsatellite  point.   The  natural  unit  of  error  meas- 
ure for  the  software,  satellite,  and  SDB  is  the  full  resolution 
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angular  interval  or  pixel.   All  elements  of  the  earth  location 
process  involving  computations  and  data  handling  have  been  de- 
signed and  implemented  to  achieve  accuracies  of  a  small  frac- 
tion of  a  pixel.   However,  the  end  to  end  process  has  several 
potential  error  sources  that  must  be  dealt  with  and  which  limit 
the  ultimately  achieveable  accuracy. 

The  potential  error  sources  are  exhaustively  listed  below. 

1.  Tolerances  in  the  orbit  and  attitude  modeling. 

2.  Tolerances  in  the  geometry  transformation. 

3.  Measurements  of  landmark  and  star  positions. 

4.  Misassignment  of  the  earth  coordinates  of  the 
landmarks. 

5.  Tolerances  on  the  satellite's  sun  detector. 

6.  Biases  in  the  SDB  scan-line  interpolation  and  sun 
position  referencing. 

7.  Timing  jitter  in  the  satellite  and  SDB. 

8.  Tolerances  in  the  solution  for  the  orbit,  attitude  and 
misalignment  parameters. 

The  timing  jitter  of  the  satellite  and  SDB  and  the  satel- 
lite's sun  detection  error  might  be  taken  as  the  hard  limits  on 
accuracy.   There  are,  however,  22  time-counts  per  full  resolu- 
tion pixel.   The  jitter  amounts  to  no  more  than  a  few  counts 
(i.e.,  less  than  a  quarter  of  a  pixel) .   The  sun  detection  er- 
ror is  the  combination  of  random  errors  in  the  leading-edge- 
detection,  which  should  be  less  than  two  or  three  counts,  and  a 
slowly  changing  offset  due  to  changes  in  the  sun's  relative 
position  in  the  sun  sensor's  slot-shaped  field  of  view.   The 
slowly  changing  offset  can  be  determined  frequently  enough  to 
be  accounted  for  in  the  "roll"  parameter  to  within  a  tolerance 
of  a  fraction  of  a  pixel.   Errors  in  the  SDB ' s  scan  line  inter- 
polation and  sun  position  reference  combined,  appear  to  be  less 
than  a  half  pixel. 

The  accuracy  of  the  geometry  formulation  is  known  to  be  well 
within  1  pixel  and  is  intended  to  be  within  0.1  of  a  pixel  but 
cannot  be  truly  verified  to  that  level.   The  orbit  model  and 
the  attitude  model  are  intended  to  be  within  1  pixel.   The 
operation  of  these  models  has  been  verified  at  SSEC  to  be 
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within  a  tolerance  approaching  1  pixel.   An  effort  has  been 
made  to  keep  these  models  commensurate  with  1  pixel  accuracy. 
Nutation  of  the  satellite's  spin  axis  can  be  modelled  and  de- 
termined, but  has  not  as  yet  been  accounted  for  in  the  earth- 
location  parameters.   Nutation,  if  not  accounted  for,  can  give 
errors  of  two  or  three  visible  lines  (peak) ,  cycling  over  about 
8  scan  lines. 

The  potential  accuracy  of  star  and  landmark  position  measure- 
ments (i.e.,  observed  image  frame  coordinates)  is  within  1 
pixel  under  normal  conditions.   Human  error  and  variations  in 
human  judgement  on  features,  can  combine  to  reduce  the  routine- 
ly achieveable  accuracy  to  about  1  pixel  rms.   The  earth  coor- 
dinates of  landmarks  can  be  verified  to  within  a  1  pixel  toler- 
ance.  The  star  coordinates  are  available  from  standard  tables 
to  an  accuracy  of  better  than  0.2  pixel  (i.e.,  0.5  arc  seconds 
in  earth  centered  celestial  coordinates) .   The  solution  for 
orbit,  attitude  and  misalignment  parameters  is  susceptible  only 
to  biases  in  the  measurements  since  it  incorporates  least- 
square-error  regression. 

When  the  errors  discussed  above  are  included,  the  resulting 
accuracy  for  a  data  user  with  a  standard  geometry  transform 
routine  that  uses  the  documented  parameters  should  be  1  pixel 
rms.   This  figure  represents  errors  over  the  entire  earth  disc 
and  throughout  the  day.   The  effect  of  spin  axis  nutation  con- 
tributes in  the  worst  of  cases  to  more  than  a  1  visible  pixel 
rms  error  in  line.   It's  contribution  to  overall  accuracy  is 
however  greatly  reduced  because  of  its  high  frequency  cycle. 
Also  the  nutation  dynamically  deforms  earth  features  but  a 
human  or  computerized  recognition  process  filters  out  much  of 
its  effect. 
2.   Actual  Results 

Figure  14  presents  the  results  of  a  determination  of  orbit, 
attitude,  and  misalignment  parameters  from  star  and  landmark 
measurements.   The  results  shown  are  not  the  best  obtainable 
and  in  fact  are  on  the  verge  of  being  unsatisfactory.   This  set 
was  selected  for  discussion  here  because  the  nature  of  the 
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various  error  contributions  is  illustrated. 

The  measured  star  and  landmark  positions  are  listed  and  iden- 
tified by  satellite  number  (SS,  West  GOES  is  #24) ,  year  (YY) 
and  Julian  (DDD)  in  the  first  column,  SSYYDDD.   The  GMT  of  the 
start  of  the  image  frames  in  which  the  particular  measurement 
where  made  is  given  as  hours-minutes-seconds  in  the  second 
column,  HHMMSS.   The  LCODE  given  in  the  third  column  identifies 
the  measurements  as; 

0  a  landmark,  in  Baja,  CA 

301  a  landmark,  in  islands  of  the  South  Pacific 

351  visible  star  #1  (|3-Orion,  "Rigel") 

353  visible  star  #3  (o'-Canis  Minor  A,  "Procyon") 

354  visible  star  #4  (a-Aquilae,  "Altair") 

355  visible  star  #5  (Y-Orion,  "Bellatr ix" ) 

The  star  a-Orion  ( "Betelgeuse" )  is  also  usually  visible  in  the 
images. 

The  fourth  and  fifth  columns  are  respectively  the  line  and 
element  residuals  (differences)  between  the  observed  image 
positions  and  the  positions  specified  by  the  solution.   The 
units  for  the  residuals  are  the  full  resolution  visible  pixel. 
The  last  two  columns  give  the  latitude  and  longitude,  in 
degrees-minutes-seconds  (DDDMMSS) ,  of  the  satellite's  subpoint. 

The  solution  was  obtained  from  the  measurements  on  day  268 
and  predicted  to  day  269.   The  prediction  is  a  linear  extrapo- 
lation in  time  of  the  attitude  and  an  evaluation  of  the  orbit 
model  with  epoch  parameters  (for  0  GMT  on  day  268)  at  the  par- 
ticular times  on  day  269.   The  solution  residuals  are  roughly 
0.5  rms  in  line  and  1.5  rms  in  element.   The  prediction  residu- 
als are  roughly  2.5  rms  in  line  and  1.5  rms  in  element.   The 
measurements  are  spread  over  nearly  24  hours  and  geographically 
dispersed  with  landmark  0  being  far  above  and  to  the  right  of 
the  subpoint  and  landmark  301  slightly  below  and  to  the  left  of 
the  subpoint.   The  different  stars  vary  in  right  ascension  by 
about  220°  and  in  declination  between  ±8°. 

The  largest  component  of  the  prediction  error,  the  "bulge"  in 
line  residuals  during  the  period  17  to  22  hours  GMT,  suggests 
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an  imprecise  solution  for  the  precesion  rate  of  the  satellite's 
attitude.   The  combination  of  all  the  other  error  sources  dis- 
cussed above  is  seen  as  the  1-  to  2-pixel  variation  that  is 
somewhat  random. 
3.   A  Chronology  of  Results  with  VIRGS 

When  the  VIRGS  became  operational  in  May,  1979,  the  orbit/at- 
titude determination  software  was  at  an  intermediate  stage  of 
its  development.   The  level  of  accuracy  achieved  in  the  opera- 
tion with  this  version  of  the  software  during  the  first  several 
months  of  operation  was  sporadic,  ranging  from  3  pixels  rms  on 
some  days  to  10  pixels  on  other  days.   On  an  average  of  once 
per  week  during  this  "break-in"  period,  severely  degraded  per- 
formance occurred  as  a  result  of  procedural  and  training  defi- 
ciences  in  dealing  with  non-routine  matters  such  as  satellite 
maneuvers. 

The  intermediate  version  of  the  software,  while  still  being 
used  in  the  operation,  was  unable  to  adequately  distinguish  the 
orbit  plane  from  the  attitude  of  the  spin  axis  using  landmark 
and  earth  edge  measurements.   By  August  1,  1979  the  software 
had  been  revised  by  Dennis  Phillips  for  the  "star  technique" 
described  in  this  document.   Preliminary  tests  of  the  revised 
software  verified  that  the  attitude  and  misalignment  parameters 
could  be  accurately  determined  from  the  star  positions  alone, 
leaving  the  landmark  positions  for  an  independent  orbit 
determination. 

During  the  months  of  August,  September,  and  October  1979,  the 
star  technique's  operational  viability  was  tested  by  Larry  Ham- 
brick.   While  having  to  work  around  the  regular  operation  on 
VIRGS,  image  sectors  containing  stars  were  ingested  from  both 
the  East  and  the  West  GOES.   The  star  positions  were  measured 
and  entered  on  the  navigation  files.   The  attitude  and  misa- 
lignment parameters  were  determined  from  the  star  measure- 
ments.  With  these  parameters  and  the  landmark  measurements 
(i.e.,  from  the  ongoing  operation),  the  orbit  was  determined. 
The  accuracy  of  the  predictions  using  these  solutions  were  then 
checked.   About  fifty  daily  navigations  were  performed  in  this 
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manner.   A  formal  demonstration  was  held  in  September,  1979. 

The  following  conclusions  are  drawn  from  the  results  of  those 
tests: 

i.  A  sufficient  set  of  stars,  for  the  star  technique,  are 

routinely  available  in  the  images, 
ii.  Using  the  star  technique,  VIRGS  has  a  stand-alone 

orbit/attitude  determination  capability, 
iii.  Several  deficiencies  remain  in  the  revised  software 

(see  next  paragraph);  however,  accuracies  better  than 

3  pixels  rms  are  routinely  achievable  on  the 

predictions,  using  the  star  technique. 
At  the  time  of  these  tests  a  discrepancy  in  "beta",  apparent- 
ly between  VIRGS  and  the  SDB,  existed  at  the  level  of  about  2 
pixels.  Also  there  was  a  discrepancy,  internal  to  VIRGS,  be- 
tween the  attitude  and  orbit  routines  as  manifested  in  their 
residual  outputs.   The  computation  of  the  attitude's  precession 
rates  was  done  manually,  a  cumbersome  task  to  do  precisely. 
Work  by  Dennis  Phillips  subsequent  to  the  tests  is  expected  to 
have  remedied  all  of  these  problems. 

A  separate  type  of  deficiency  that  had  to  be  dealt  with  dur- 
ing these  tests  was  a  discrepancy  in  the  earth  coordinates  of 
landmarks.   Landmarks  not  on  the  American  Continents  appeared 
to  have  inconsistent  earth  coordinates  with  respect  to  American 
landmarks.   The  coordinates  for  all  land  features  had  been  ob- 
tained from  official  detailed  maps.   An  extensive  test  was  con- 
ducted by  Larry  Hambrick  on  VIRGS  (using  the  star  technique)  to 
establish  the  size  of  the  discrepancy  for  about  a  dozen  non- 
American  land  features.   For  some  islands  in  the  South  Pacific, 
the  discrepancy  was  as  high  as  five  image  pixels.   A 
discrepancy  of  at  least  two  pixels  was  found  for  most  of  the 
non-American  features.   Based  on  this  test  it  was  concluded 
that  a  calibration   of  the  coordinates  of  non-American 
features  with  respect  to  American  land  features  was  needed  and 
would  further  improve  the  accuracy  of  VIRGS. 
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V.   SUMMARY 

The  accurate  earth  location  of  GOES/VISSR  image  data  requires 
a  complete  mathematical  formulation  that  is  founded  on  the 
basic  geometry.   The  use  of  stars  together  with  landmarks  as 
the  observables  in  the  images  affords  a  small,  efficient  soft- 
ware package  that  employs  relatively  simple  algorithms  and 
models  for  determining  the  true  misalignment  angles,  attitude, 
and  orbit. 

The  capability  described  in  this  document  has  been  implemen- 
ted on  VTRGS,  the  central  VISSR  earth-location  facility  oper- 
ated by  NOAA/NESS  and  on  the  Mcldas  at  the  University  of  Wis- 
consin.  Direct  readout  users  have  convenient  access  to  the 
VIRGS  generated  earth-location  parameters,  and  by  using  a  stan- 
dard geometry  transform  routine,  can  locate  the  VISSR  data 
automatically.   The  NESS  standard  grids,  also  produced  by 
VIRGS,  are  embedded  in  the  VISSR  data  at  the  infrared  pixel 
resolution. 

The  VIRGS  should  be  operated  at  the  best  achievable  level  of 
accuracy.   A  commitment  to  an  accurate,  reliable  centralized 
service  would  save  individual  users  a  sizable  task.   The  users 
face  the  choice  of  exercising  their  own  VIRGS-type  operation  or 
simply  using  the  earth-location  parameters  of  the  centralized 
service. 
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APPENDIX  A 
CONVERTING  A  POINTING  VECTOR  TO  LINE  AND  ELEMENT  NUMBERS 
The  pointing  vector  in  the  first  satellite  centered  coordi- 
nate system  is  designated  as  k.   The  projection  of  k  onto  the 
z-axis  of  the  first  satellite  coordinate  system  is  simply 
k»(-S),  where  S  is  the  spin  axis  vector.   In  terms  of  scan 
line  number  and  the  misalignment  parameters,  the  projection  of 
the  pointing  direction  of  the  spin  scan  camera  onto  the  z-axis 
of  the  first  satellite  coordinate  system  is  simply  the  z  com- 
ponent of  (2.5)  or 

m^  ,  ocos  (  (c*  -  H)  .  r,  )  +  m3  3#sin(  (Cn-{)  t»)  .       (A.O) 
To  determine  the  line  number  at  which  the  spin  scan  camera  is 
pointing  in  the  direction  k,  simply  equate  (A.O)  to  k»(-S)  and 
solve  for  12  .   Solving  this  equation  requires  a  trigometric 
manipulation. 

2      2  Vi 
Divide  (A.O)  by  (m3  1   +ir.3  3  )  , 

set  9=ATAN2(m3  x,m3  3)  and 

2      2/2 
substitute  sin(0)  for  itu  ,/(m3  .  +mr.    3  )   and 

2  '      2  !/2 
cos  (9)  for  m3  3/(m3  ,  +m->  3  )  . 

The  result  is 

sin  (9)  o  cos(  (c^,  -St )  «r«  )  +  cos  (9)  •  sin(  (c»  - H)  •  r*  ) 

=r-  (-"S)/dn3fl+m3^3)1/2 

Setting  4>  =  ASIN  (k.  (-"§)  /  (m^  x+m3  3)  2 )  , 
using  the  sum  angle  sine  angle  formula  on  the  left-hand  side  and 
applying  the  arcsine  function  to  both  sides  gives: 

(Cg-2)  .r^  +  0  =  <K  (A.l) 

This  linear  expression  for  line  number  is  then  solved.   The 
resultant  value  of  the  line  number  is  checked  to  ensure  it  lies 
within  the  bounds  for  the  image,  i.e.,  1^£<14568. 

Finally,  the  element  count  is  to  be  found.   This  is  done  by 
equating  the  angular  positions  in  the  spin  plane,  around  the 
z-axis  of  our  first  satellite  coordinate  system,  as  given  by 
the  camera  geometry  and  as  given  by  the  k  projection.   By 
taking  arctangents  of  the  ratio  of  the  x  and  y  components  of 
the  vector  formulated  by  (2.5)  one  finds  that  the  angular 
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component  of  this  vector  in  the  spin  plane  is  given  as  follows: 

-(Y-(3.tp.sr  -  s.rg)  + 
ATAN2[m2  ,cos(\)+m2  3sin(\):m,  ,cos(\)+m,  3sin(\)]  (A. 2) 

where  X.  =  {  (c»-H)  ,rn)  • 

The  x  component  and  y  component  of  the  pointing  vector  k  in 
the  same  coordinate  system  are  found  respectively  by  projec- 
tions k»P  and  k»PXS_  where  P  and  S_  come  from  the  matrix 

P  P 

formulated  in  (2.6).   Hence,  the  angular  position  of  this 
pointing  vector  in  the  spin  plane  equals  ATAN2(k»PXS  ,  k»P) . 

IT 

This  is  simply  equated  to  (A. 2)  and  the  resultant  linear 
equation  is  solved  for  the  sample  number.   Some  care  is 
required  so  that  an  incorrect  multiple  of  2  tt  is  not 
inadvertently  included  in  one's  answer. 
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APPENDIX  B 
DETERMINATION  OF  ATTITUDE  AND  ORBITAL  PLAN  ORIENTATION 
In  order  to  determine  the  satellite's  attitude  and  the  ori- 
entation of  the  satellite's  orbital  plane  one  has  the  task  to 
minimize  the  following  expression: 

k  2 

S(a,  b,  c)=  S  (ax  .+by  .+cz  .  -cos(  cj>.)  )  (B.l) 

i=l 

2   2   2 
subject  to  the  constraint  a  +b  +c  =1  and  where  for  atti- 
tude determination  cj>  .  equals  TT/2+r  -.  ( c,  -1  . )  ,  see  figure 
8,  and  for  the  orientation  of  the  orbital  plane  all  the  <j>  ,  '  s 

equal  1T/2.   The  constraint  is  absorbed  by  setting  c  = 

2   2  Vo 
-(l-(a  +b  )  )n  .   Starting  with  a  and  b  equal  to  zero  sets 

the  initial  spin  axis  vector  guess  and  the  orbital  plane  per- 
pendicular guess  pointing  opposite  to  the  spin  axis  of  the 
earth.   This  is  a  good  operational  approximation  and,  conse- 
quently, reduces  the  number  of  iterations  required  during 
computation.   The  change  of  variables  modifies  the  problem  to 
minimizing  the  following: 

S(a,b)  =2  (ax.+byi-(l-(a2+b2)  )/2.zi-cos(4>i)  )  2.     (B.2) 
i=l 
The  first  partials  of  S  with  respect  to  a  and  b  are  set  equal 

to  zero  resulting  in  the  following  two  equations: 

k 

2  (x^a^  aJ  • (axi+byi+czi-cos(  4^) ) =0 

1^1  (B.3) 

2  (y.-bz./  ) . (ax . +by . +cz . -cost §• ) ) =0 
._,  J  1    l/c       1  2  1    1        1 

2      2     ]A 

Note    the   use   of    the    symbol   c    here    for    -(l-(a    +b    ))     .      Ma- 
nipulation  yields    the    following   matrix   equation   for    the    indi- 
vidual   terms   of    the    summation: 

2        2 
x.    -z.    +z  .cos(  <)>.  ) /c  x.y. 

1111  iJ  l 


2         2  ,       cos  (<b  . ) 
x.y.  y.    -z.    +z yi 

i-1  l  2  l         l         l        r 


(B.4) 


2 

x.cos(<t>.)-cx.z.+a   x.z.  /c+aby  .  z  .  /c 
l  l  11  11  2  i    l 

2 

y.cos(4,-)-cy.z.  +abx  .  z  .  /c+b   y  .  z  .  /c 

jr  1  t1/         :L     1  x     ■}_'  J  i     i' 
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The  matrix  on  the  left  is  inverted  and  the  resulting  system  of 
equations  is  solved  iteratively.   Note  that  all  the  terms  on 
the  right-hand  side  outside  of  the  constant  terms  are  second 
order.   The  advantage  of  such  a  separation  is  that  the  linear 
equation  is  solved  for  first  and  then  the  results  are  updated 
with  the  second  order  terms.   This  approach  is  a  variation  of 
Newton's  Method.   The  summation  signs  have  been  omitted  to  com- 
pact the  presentation.   Convergence  is  achieved  within  three  or 
four  iterations. 

For  attitude  determination  with  an  unknown  pitch  misalignment 
angle, C,  the  following  expression  of  three  variables  has  to  be 
minimized: 

k  2 

S(a,  b,  c,  £  )=  s  (ax  +by  +cz  -cos(4>  -£))         (B.5) 

,111      i 
i=l 

2   2   2 
subject  to  the  constraint  a  +b  +c  =1. 

The  expression  (B.5)  is  not  easily  amenable  to  the  variation 

of  Newton's  Method  used  earlier  to  find  the  attitude  and  hence 

the  expression  is  reformulated  to  an  equivalent  problem.   The 

*2   2   2 
constraint  a  +b  +c  =1  is  absorbed  by  setting 

2   2  xh  2  !/2 

c=  -(l-(a  +b  ))   and  u  is  set  to  sin(C)  and  -(1-u  )   is 

set  to  cos(C) •   Once  u  is  solved  for,  C  can  be  found  as  the 

2  V2 
arctangent  of  (-u/(l-u  )  )  ,  or  in  Fortran  as 

2  Vi 
ATAN2 (u,- (1-u  )')•   With  this  formulation,  the  following  is 

minimized: 

S(a,b,u)=  Z  (ax.+by.-(l-(a2+b2)  z.+ucos(4>.)  +  (l-uV2sin(4>.  )  )  2. 

1=1  (B-6) 

To  solve  this  minimization  problem  the  first  partials  of  S  with 
respect  to  a,  b,  and  u  are  set  equal  to  zero  and  the  resulting 
expressions  are  rearranged  to  have  a  matrix  times  the  vector 
(a,  b,  u)  on  the  left-hand  side  and  constant  terms  plus  second 
order  terms  on  the  right-hand  side",  as  in  (B-4)  .   The  inverse 
of  the  matrix  is  then  applied  to  both  sides  and  the  resulting 
equations  are  solved  by  iteration  using  Newton's  Method.   Con- 
vergence is  again  achieved  within  three  or  four  iterations. 
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APPENDIX  C 
ACCURACY  OF  APPROXIMATE  ORBIT  EQUATION 
(Refer  to  Section  III-D,  Equation  4.4) 
In  this  paper  the  equation 

n(t-T)=  Y  -  2esinY  (C.l) 

has  been  used  to  approximate  the  motion  of  a  satellite  in  a 
Keplerian  orbit.   Equation  C.l  is  deriveable  from  (4.4)  with 
the  substitution   Y  =  a   -  co.   In  eq.  (C.l)  ,  n  is  the  mean 
motion  constant,  t  is  time,  T  is  the  time  of  a  perifocus 
passage,  Y  is  the  true  anamoly  and  e  the  eccentricity.   The 
accuracy  of  this  equation  for  eccentricities  less  than  or  equal 
to  0.001  will  now  be  established. 

Consider  the  following  three  equations  taken  from  Escabol: 

n(t-T)  =  E  -  €  sin  E,  (C.2) 


=  A--2 


sin  E   =  vi-€  sinY  (C.3) 

1+fcosY         and 

cos  E   =  cos  Y  +  £  (C.4) 

1+ecos  Y 

where  E  is  the  eccentric  anamoly. 

Multiplying  (C.3)  by  cos Y  and  (C.4)  by  sin Y and  subtracting 

yields;   cos Y  sin  E  -  sin Y cos  E  = 
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e  sinY  cos  Y  t  sinYcos  Y  -  e  »  sinY  (C.5) 

1  +  c • cos  Y 

Using  the  sine  angle  difference  formula  changes  that  left  hand 

side  to  sin  (E-Y) . 

The  next  two  inequalities  are  well  know  for  |E-Y|  within  the 

range  of  consideration: 

Isin  (E-Y)|>|E-Y|  (l-(E-Y)2) 

3!  and         (C.6) 

Isin  (E-Y)  -  (E-Y) 1<1E-Y) 3) 

3!  (C.7) 

Using  (C.5)  and  restricting  e  to  be  less  than  or  equal  to  0.001 
yields 

Isin(E-Y) |<0. 001002. 
Combining  this  with  (C.6)  yields 

|E-Y|<. 001003 
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Finally,  (C.7)  yields 
I  sin  (E-Y)  -  ( 
Now  (C5)  is  rewritten  as 


Isin  (E-Y)  -  (E-Y) |< 1. 682xl0~10  (C.8) 


r~  2  2 

vl-e    sinYcosY-sinYcosY+  €     sinYcosY 

sin(E-Y)=-«sinY+ =— r ^7— — — 

1    +  f»cos  Y 

Combining  (C.9)  with  the  inequality  (C.8)  and  appropriately 


bounding  the  last  term  in  (C.9)  yields 
|E-Y  +  e  sinYl  <  1.6  x 
Next,  (C.3)  is  rewritten  as 


+  e  sinYl  <  1.6  x  10"6  (C„ 10) 


I        2 

vl-  e   sin  V-  sin  Y  - e  »  cosYsinY       /r.  ,  -,  v 

„  •   v \  C  •  1 1 ) 

sin  E  -  sin  «  =  r — tt~~ -—-—_ — _____ 

1  +  e  .  cos  Y 

This  yields  the  inequality 

Isin  E  -  sinYl  <  .001002  (C.12) 

Combining  (1.2)  with  the  inequalities  (CIO)  and  (C.12)  yields 
In(t-T)  -  Y  +  2e  sinYl  <  2.7  x  10~6  (C.13) 

This  yields  a  positional  accuracy  to  within  0.12  kilometers. 
Thus  the  accuracy  of  the  approximation  of  using  (C.l)  is 
established. 
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APPENDIX  D 

I  GAUSSIAN  ELIMINATION 

The  process  of  Guassian  elimination  is  now  recalled  for  the 

reader.   One  begins  with  a  matrix  pair  (A:I) ,  where  A  is  the 

matrix  whose  inverse  is  sought  and  I  is  the  identity  matrix, 

and  applies  a  sequence  of  non-singular  matrices   B.   to  both 

members  of  the  matrix  pair  from  the  left  until 

B  B   ,...B,A  is  the  identity  matrix.   At  that  point  the 
n  n-1     1  ■* 

product  accumulated  on  the  right,  B  B   n...B, ,  is  the 

n  n-1     1 

inverse  of  the  matrix  A  since  its  product  with  A  yields  the 
identify  matrix.   Usually  the  B.'s  are  selected  so  that  the 
matrix  multiplication  affects  only  one  row  of  the  matrix  pair 
Because  the  matrix  (4.12)  is  symmetric,  positive  definite, 
i.e.,  a.  .=a.  .,  and  all  the  eigenvalues  are  positive,  the 
matrix  operands  B.   can  be  selected  to  act  on  more  than  one 
row  in  one  operation.   In  this  case  the  B.'s  act  on  two 
rows.   The  matrix  (4.12)  is  treated  as  if  it  were  partitioned 
into  four  parts,  each  part  consisting  of  a  two  by  two  matrix. 
Gaussian  elimination  is  then  applied  as  if  each  two  by  two 
submatrix  was  a  single  element. 


Start  with  the  matrix  pair  /A,  ,  A,  0\     /I   0\  and  apply 


Next  apply/  I   -A,  2  \  with  the  result 


Vo 


I 


Al,l"Al,2A2,2   A2,l    °\      I1         ~A1,2A2,2 


A2,2   A2,l        V   \°      A2,2 


D-l 


Then  set  S  =  (A,  ,  -  A1  2A2  „  A^      )         and  apply  /S   O 

_1\o     I 
A1,2A2,2 

"I 

A2,2 


Finally >    apply/      I     O\to  both  matrices  with  the  result 

-1 
"A2,2   A2,l 


S  "SA1,2A2,2 

-1  -1  -1  -1 

~A2,2      A2,1S  A2,2      +A2,2      A2 , 1SA1, 2A2 , 2 

The  inverse  of  (4.12)  is  the  right-hand  term.   Actually,  com- 
puting the  lower  left-hand  term  is  unnecessary  since  it  is  the 
transpose  of  the  upper  right  term  by  symmetry. 
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Appendix  E 

Fortran  Listings  of  Orbit  and  Attitude 

Determination  Routines: 

FINDAT  and  UPGORB 
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SUBROUTINE  FINDAT( NOPCLN .RADLIN .DECLIN .RASCEN , PCLN ) 

DOUBLE  PRECISION  UX.UY.UZ.PSI 

COMMON/MINCOM/UX(60) ,UY (60 ) , UZ (60 ) , PS  I ( 80) , NP 

DATA  RDPDG/1.745329252E-2/ 
C     WRITTEN  NOVEMBER  13,1979  BY  DENNIS  PHILLIPS  TO  SIMPLIFY  ATTITUDE 
C     COMPUTATION.   PROGRAM  RIGHTS  BELONG  TO  SCIENTIFIC  PROGRAMMING 
C     AND  APPLIED  MATHEMATICS  ,  INC. 

C     THIS  PROGRAM  ELIMINATES  THE  USE  OF  STEEPEST  DESCENT  METHODS  IN 
C     THE  ATTITUDE  COMPUTATIONS  AND  INSTEAD  REPLACES  THESE  METHODS  WITH 
C     LINEAR  REGRESSION  METHODS  COMBINED  WITH  ITERATION  STEPS. 
C     TWO  PROBLEMS  ARE  WORKED  OUT  IN  THIS  CODE.   ONE  PROBLEM  IS  TO 
C     MINIMIZE  THE  SUM  OVER  I  OF  (A*X ( I )  +  B*Y ( I )+C*Z ( I )-COS (PSI (I  )  )  )**2. 
C     WITH  THE  CONSTRAINT  A**2+B**2+C**2=l . 
C     THE  OTHER  PROBLEM  IS  TO  MINIMIZE  THE  SUM  OVER  I  OF 
C      (A*X(I.)+B*Y(I)  +  C*Z(I)-C0S(PSI(I)-PSIREF))**2  OR  FQUIVALENTLY  THE 
C     SUM  OVER  I  OF  ( A*X( I )+B*Y ( I ) +C*Z ( I )-V*COS ( PSI ( I ) ) -U*SIN ( PSI ( I ) ) )**2 
C     WITH  THE  CONSTRAINTS  A**2+B**2+C**2  AND  U**2+V**2=l. 
C     *** 

C     FOR  PICTURE  CENTER  LINE  SOLUTION  GO  TO  30 

Q  ### 

IF(NOPCLN.EQ.l)GO  TO  30 
C     RETURN  AND  SEND  ERROR  MESSAGE  FOR  INSUFFICIENT  MEASUREMENTS 
IF(NP.LT.2)G0  TO  75 

C     COMPUTE  ATTITUDE 
C     *#* 

C     COLLECT  REGRESSION  COEFFICIENTS 

SMXSQ=0.0 

SMXY=0.0 

SMXZ=0.0 

SMXCS=0.0 

SMYSQ=0.0 

SMYZ=0.0 

SMYCS=0.0 

SMZSQ=0.0 

SMZCS=0.0 

DO  10  1=1, NP 

X=UX(I) 

Y=UY(I) 

Z=UZ(I) 

ANG=PSI(I) 

CS-COS(ANG) 

SMXSQ=SMXSQ+X*X 

SMXY=SMXY+X*Y 

SMXZ=SMXZ+X*Z 

SMXCS=SMXCS+X*CS 

SMYSQ=SMYSQ+Y*Y 

SMYZ=SMYZ+Y*Z 

SMYCS=SMYCS+Y*CS 

SMZSQ=SMZSq+Z*Z 
10  SMZCS=SMZCS+Z*CS 
C     SET  INITIAL  VALUES  AND  UPDATE  VALUES 

A=0.0 
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B  =  0„0 

C— 1.0 

ANEW=0.0 

BNEW=0.0 

SET  ITERATION  COUNT 

ITER=0 

COMPUTE  MATRIX  COEFFICIENTS  FOR  MATRIX  EQUATION 


(All 
(A21 


A12) 
A22) 


=    Fl 
=    F2 


20   A11=SMXSQ+SMZCS/C-SMZSQ 

A12-SMXY 

A21=SMXY 

A22=SMYSQ+SMZCS/C-SMZSQ 
C      INVERT  THE  MATRIX  A 

DET=1.0/(A11*A22-A12*A21 ) 

B11=A22*DET 

B12=-A12*DET 

B21=B12 

B22=A11*DET 

COMPUTE  CONSTANT  PLUS  SECOND  ORDER  TERMS  OF  THE  VECTOR  F 

F1=SMXCS-C*SMXZ+ANEW**2*SMXZ/C+ANSW*BNEW*SMYZ/C 

F2=SMYCS-C*SMYZ+BNEW**2*SMYZ/C+ANEW*BNEW*SMXZ/C 
C  COMPUTE    NEW    A    AND    B   VALUES 

ANEW=B11*F1+B12*F2 

BNEW=B21*F1+B22*F2 

C=-SQRT(1.0-(ANEW**2+BNEW**2 ) ) 
C      COUNT  ITERATION  STEPS 

ITER=ITER+1 
C     RETURN  AND  SEND  ERROR  MESSAGE  IF  ITERATION  LIMIT  IS  EXCEEDED 

IF(ITER.GE.10)GO  TO  80 
C      CHECK  CONVERGENCE  CRITERIA 

IF(ABS (A-ANEW)+ABS(B-BNEW) .LT.1.0E-8)GO  TO  25 
C      UPDATE  VALUES 

A=ANEW 

3=BNEW 
C      ITERATE  AGAIN 

GO  TO  20 
25  DECLIN=-90.0+ASIN(SQRT(A**2+B**2))/RDPDG 

RASCEN=0.0 

IF(ABS(A)+ABS(B) .GT . 1 . 0E-8 )RASCEN=ATAN2 (B , A) /RDPDG 
RETURN 
C     *#* 

C      COMPUTE  ATTITUDE  AND  PICTURE  CENTER  LINE  OFFSET 
C      *** 

30  IF(NP.LT.3)G0  TO  85 
C      COLLECT  REGRESSION  COEFFICIENTS 
SMXSQ-0.0 
SMXY=0.0 
SMXZ=0.0 
SMXCS-0.0 
SMXSN=0.0 
SMYSQ=0.0 
SMYZ=0.0 
SMYCS-0.0 
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SMYSN=0.0 

SMZSQ-0.0 

SMZCS=0.0 

SMZSN=0.0 

SMCSSQ=0.0 

SMCSSN=0.0 

SMSNSQ=0.0 

DO  35  1=1, NP 

X  =  UX(I)' 

Y=UY(I) 

Z=UZ(I) 

ANG=PSI (I) 

SN=SIN(ANG) 

CS-COS(ANG) 

SMXSQ=SMXSQ+X*X 

SMXY=SMXY+X*Y 

SMXZ=SMXZ+X*Z 

SMXCS=SMXCS+X*CS 

SMXSN=SMXSN+X*SN 

SMYSQ=SMYSQ+Y*Y 

SMYZ=SMYZ+Y*Z 

SMYCS=SMYCS+Y*CS 

SMYSN-SMYSN+Y*SN 

SMZ3Q=SMZSQ+Z*Z 

SMZCS=SMZCS+Z*CS 

SMZSN=SMZSN+Z*SN 

SMCSSQ=SMCSSQ+CS*CS 

SMCSSN=SMCSSN+CS*SN 
35  SMSNSQ=SMSNSQ+SN*SN 
C      SET  INITIAL  AND  UPDATE  VALUES 

A  =  0.0 

B=0.0 

C=-1.0 

ANEW=0.0 

BNEW=0.0 

U=0.0 

V  =  1.0 

UNEW=0.0 
C     SET  ITERATION  COUNT 

ITER=0 
C     COMPUTE  COEFFICIENTS  FOR  MATRIX  EQUATION  (All  A12  A13)  A  =  Fl 
C  (A21  A22  A23)  B  =  F2 

C  (A31  A32  A33)  U  =  F3 

40  A11=SMXSQ-SMZSQ+V*SMZCS/C 

A12=SMXY 

A13=-SMXSN 

A21-A12 

A22=SMYSQ-SMZSQ+V*SMZCS/C 
t    A23=-SMYSN 

A31-A13 

A32=A23 

A33=C*SMZCS/V-SMCSSQ+SMSNSQ 
r  #** 
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C  THE  INVERSION  OF  THE  MATRIX  A  IS  DONE  BY  PARTITIONING  THE  MATRIX 

C  INTO  FOUR  PARTS:  (All),  (A12  A13)  ,  (A21)  AND  (A22  A23) 

C  (A31)      (A32  A33) 

C  THEN  EACH  MEMBER  OF  THE  PARTITION  IS  TREATED  AS  A  SINGLE  ELEMENT 

C  AND  GAUSSIAN  ELIMINATION  IS  CARRIED  OUT  STARTING  IN  THE  LOWER  LEFT 

C  HAND  CORNER 

Q  *** 

C      FIND  THE  INVERSE  OF  THE  LOWER  LEFT  HAND  CORNEP 

DET=1.0/(A22*A33-A23*A32) 

CU=A33*DET 

C22=A22*DET 

C12=-A23*DET 

C21=C12 
C      MULTIPLY-  THIS  INVERSE  BY  THE  LAST  TWO  ENTRIES  OF  THE  FIRST  ROW 

B1=-A12*C11-A13*C21 

B2=-A12*C12-A13*C22 
C     FIND  FIRST  ELEMENT  OF  INVERSE  MATRIX 

CDIVID=1.0/(A11+B1*A21+B2*A31 ) 

B11=CDIVID 
C      FIND  THE  SECOND  AND  THIRD  ELEMENTS  OF  FIRST  ROW  OF  MATRIX 

B12=B1*CDIVID 

B13=B2*CDIVID 
C      SET  TWO  MORE  ENTRIES  BY  SYMMETRY 

B21=B12 

B31=B13 
C      COMPUTE  THE  REMAINING  TWO  BY  TWO  PART 

B22=C11+CDIVID*B1*B1 

B23=C12+CDIVID*B1*B2 
C     BY  SYMMETRY 

B32=B23 

B33=C22+CDIVID*B2*B2 
C      COMPUTE  RIGHT-HAND  SIDE  OF  LINEAR  EQUATION 

F1=V*SMXCS-C*SMXZ+A**2*SMXZ/C+A*B*SMYZ/C-A*U*SMZSN/C 

F2=V*SMYCS-C*SMTZ+A*B*SMXZ/C+B**2*SMYZ/C-B*U*SMZSN/C 

F3=-V*SMCSSN+C*SMZSN-A*U*SMXCS/V-B*U*SMYCS/V+U**2*SMCSSN/V 
C      COMPUTE  NEW  A,  B  AND  U  VALUES 

ANEW=B11*F1+B12*F2+B13*F3 

BNEW=B21*F1+B22*F2+B23*F3 

UNEW=B31*F1+B32*F2+B33*F3 

ITER=ITER+1 
C      RETURN  AND  SEND  ERROR  MESSAGE  IF  ITERATION  LIMIT  IS  EXCEEDED 

IF(ITER.GE.15)G0  TO  90 
C      CHECK  CONVERGENCE  CRITERIA 

IF(ABS(A-ANEW)+ABS(B-BNEW)+ABS(U-UNEW) .LT . 1.0E-8)GO  TO  45 
C     UPDATE  VALUES 

A=ANEW 

B=BNEW 

C=-SQRT(1.0-(A*A+B*B)) 

U=UNEW 

V=SQRT(1.0-U**2) 
C      ITERATE  AGAIN 

GO  TO  40 
C      COMPUTE  DECLINATION,  RIGHT  ASCENSION  AND  PICTURE  CENTER  LINE  OFFST 
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45  DECLIN=-90.0+ASIN(SQRT(A*A+B*B)  )/RDPDG 

RASCEN=0.0 

IF(ABS(A)+ABS(B).GT.1.0E-8)RASCEN=ATAN2(B,A)/RDPDG 

PCLN=PCLN-ATAN2(U,V)/RADLIN 

RETURN 
C      SEND  ERROR  MESSAGES 
75  WRITE(6,95) 

95  FORMAT (5X, 'INSUFFICIENT  MEASUREMENTS  FOR  ATTITUDE  DETERMINATION') 
RETURN 

80  WRITE(6,96) 

96  F0RMAT(5X, 'ITERATION    LIMIT   EXCEEDED    IN    ATTITUDE   DETERMINATION') 
RETURN 

85  WRITE(6,97) 

97  F0RMAT(5X, 'INSUFFICIENT  MEASUREMENTS  FOR  ATTITUDE  AND  PCLN  DET . ' ) 
RETURN 

90  WRITE(6,98) 

98  F0RMAT(5X, 'ITERATION  LIMIT  EXCEEDED  IN  ATTITUDE  AND  PCLN  DET.') 
RETURN 

END 

SIZE     945    01661 
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C  UPGORB  PHILLI  067S  NAVL  COMPUTE  ORBIT  FROM  LANDMARKS , BETAS ,  AND  ATTITU0004 

C 

C      THIS  PROGRAM  HAS  TWO  MODES:  ONE  MODE  TO  COMPUTE  ALL  THE  ORBIT 

C      PARAMETERS  AND  THE  OTHER  MODE  TO  COMPUTE  JUST  THE  ALONG-TRACK 

C      ORBIT  PARAMETERS  ASSUMING  THE  ORBIT'S  INCLINATION  AND 

C  ASCENDING  NODE  HAVE  BEEN  PROVIDED.   MIN(3)  CONTROLS  THE 

C  NUMBER  OF  ALONG-TRACK  ORBIT  PARAMETERS  TO  BE  COMPUTED. 

C      MIN(6)  IS  THE  FLAG  CONTROLING  THE  ORBITAL  PLANE  CALCULATION. 

C 

C      MIN(l)  =  SSYYDDD 

C      MIN(2)  =  HHHHt  THE  FIRST  TWO  DIGITS  ARE  THE  FIRST  HOUR  ON  THE  YYDD 

C  TO  START.   THE  LAST  TWO  DIGITS  ARE  THE  NUMBER  OF 

C  HOURS  FROM  THE  FIRST  HOUR  TO  CONSIDER.  0 

C      MIN(3)  =  NO.  OF  PARAMETERS  TO  COMPUTE. 

C      MIN(4)  =  EPOCH  OF  THE  CALCULATED  ORBIT  PARAMETERS.  00 

C      MIN(5)  =  PRINT  AND  STORE  OPTION. 

C      MIN(6)  =  ORBITAL  PLANE  CALCULATION  FLAG.   OP  MEANS  ON. 

SUBROUTINE  MAIN  0012 

REAL  MEANOM,PTIMLM(100) ,XLINLM(100) fXELELM(100)  0013 

REAL  XLATLM(100) ,XLONLM(100)  0014 

DOUBLE  PRECISION  BETA( 100 ) ,V ( 100) ,T (100 )  0015 

INTEGER  HHHH  0016 

INTEGER  MIN(S) ,LCODE(100) ,ILNDBT(100)  0017 

INTEGER  LNDPNT(100)  0018 

COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET0019 
1IMH,SEMIMA,OECCEN,ORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN0020 
2, PRERAT.PREDIR, PITCH, YAW, ROLL, SKEW  0021 

DATA  MIN/'UPGORB',6*0/  0022 

DATA  MAXLND/100/,NUMDAY/1/  0023 

CALL  IQ(MIN)  0024 

ISYD=MIN(1)  0025 

HHHH=MIN(2)  0026 

NUMDAYMHHHH/100+MOD  (HHHH,  100)  )/24+l  0027 

C      THE  CALL  TO  SETUPN  INITIALIZES  THE  PARAMETERS  IN  THE  NAVCOM  AND 

C      NAVIN  COMMON  BLOCKS 

CALL  SETUPN(ISYD,0)  0028 

C     FETCH  THE  LANDMARKS  FOR  THE  TIME  PERIOD 

CALL  GTLM(ISYD,NUMDAY,MAXLND,NUMLND,  0029 

*  PTIMLM,LCODE,XLINLM,XELELM,XLATLM,XLONLM)  0030 
C     THIS  CALL  TO  GTLMBT  MATCHES  BETA  VALUES  TO  THE  AVAILABLE  LANDMARKS 

C      BY  INTERPOLATING  OR  EXTRAPOLATING  THE  BETA  VALUES  THAT  ARE 

C     AVAILABLE  IN  THE  SAME  IMAGE  WITH  THE  LANDMARKS.   WHEN  SUCH  BETAS 

C     ARE  AVAILABLE  FOP  A  SPECIFIC  LANDMARK,  THE  CORRESPONDING  INDEX 

C      POSITION  OF  THE  LANDMARK  IN  THE  ARRAY  ILNDBT  IS  SET  TO  ONE? 

C      OTHERWISE  THIS  VALUE  IS  SET  TO  A  MINUS  ONE. 

CALL  GTLMBT ( I SYD,NUMDAY, HHHH, NUMLND ,PTIMLM .XLINLM ,XELELM , LCODE, 

*  ILNDBT, BETA, T) 

C      COUNT  LANDMARKS  HAVING  ASSOCIATED  BETA  VALUES 

IMATCH=0  0033 

DO  2  1=1, NUMLND  0034 

IF( ILNDBT (I ) .EQ.1)IMATCH=IMATCH+1  0035 

2     CONTINUE  0036 

C      TRANSMIT  ERROR  MESSAGE  AND  RETURN  WHEN  NOT  ENOUGH  LANDMARKS  HAVE 
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ASSOCIATED  BETA  COUNTS 

IF(IMATCH.GE.MIN(3))G0  TO  3  0037 

CALL  TQMESCNOT  ENOUGH  MATCHES  BETWEEN  LANDMARKS  AND  BETAS $', IMATC003S 

*H)  0039 

RETURN  0040 

CONTINUE  0041 

IF(MIN(3).GE.1.AND.MIN(3) .LE.4)GO  TO  4  0042 
CALL  TQMES(  'WRONG  NUMBER  OF  ORBIT  PARAMETERS  TO  BE  COMPUTED??',     0043 

*MIN(3))  0044 

RETURN  0045 

ISEMI=MIN(3)  0046 

EPTIME=FTIME(MIN(4))  0047 

IHNS=MIN(4)  0048 

IORBPL=0  0049 

IF(MIN(6),EQ.2H0P)IQRBPL=1  0050 

THE  ITERATION  LOOP  (DO  6  1=1,3)  CORRECTS  THE  HEIGHT  OF  THE 
SATELLITE  AS  A  FUNCTION  OF  TIME  USING  THE  MOST  RECENTLY 
COMPUTED  VALUES  FOP  THE  ORBIT  PARAMETERS.   WHEN  NO  ORBIT 
PARAMETERS  ARE  AVAILABLE,  THE  NOMINAL  GEOSTATIONARY  HEIGHT  OF 
42165.0  IS  USED.   THIS  HEIGm.  CORRECTION  IS  DONE  IN  ANGORB 

IETIMH=MIN(4) 

IETIMY=MOD( I SYD, 100000) 

ITER=5 

MEANOM=0.0 

DO  6  1=1, ITER 

THE  CALL  TO  ANGORB  SERVES  TWO  FUNCTIONS:  (1)  IF  IORBPL=l ,  THE 
INCLINATION  AND  ASCENDING  NODE  OF  THE  ORBITAL  PLANE  ARE  DETERMINED 
AND  (2)  THE  TRUE  ANAMOLY  POSITIONS  +  OFFSET  OF  THE  SATELLITE  ARE 
COMPUTED.   THE  TRUE  ANAMOLIES+OFFSET  ARE  STORED  IN  THE  ARRAY  V(I). 

CALL  ANGORB ( I SYD, NUMLND ,MEANOM ,PTIMLM ,XLATLM,XLONLM , 
*XLINLM,XELELM,ILNDBT,BETA,V,T,NV,LNDPNT,IORBPL)  0059 

COMPUTE  THE  ALONG-TRACK  ORBIT  PARAMETERS.   THE  NUMBER  OF 
PARAMETERS  COMPUTED  DEPENDS  ON  THE  INDEX  ISEMI .   HENCE  I  SEMI 
EQUALS  ONE,  TWO,  THREE  OR  FOUR. 


0054 
0055 


0052 


CALL  ORBPR(V,T,NV,MEANOM, ISEMI, EPTIME) 
6  CONTINUE 
IOUT=l 

COMPUTE  RESIDUALS  FOR  LANDMARK  MEASUREMENTS  USING  THE  NEWLY 
COMPUTED  ORBITAL  STATE  AND  OUTPUT  RESIDUALS  TO  TERMINAL. 

CALL  ORBRES(ISYD,IOUT,NV,T,XLINLM,XELELM,BETA,MEANOM, 
*     LCODE,XLATLM,XLONLM,LNDPNT) 
I0UT=2 

OUTPUT  RESIDUALS  TO  PRINTER  IF  PRINTER  OPTION  IS  ON. 


0060 
0061 
0062 


0064 
0065 


18.050575  PAGE       3 

IF(MIN(5).FQ.2HPS.0R.MIN(5) .EQ . 2HSP .OR . MI N ( 5 ) .EQ. 1HP ) CALL  ORBRES   0066 
*(ISYD,IOUT,NV,T,XLINLM,XELELM,BETA,MEANOM,LCODE,XLATLM.XLONLM,LNDP 
*NT) 
C 

C      STORE  NAVIGATIONAL  PARAMETERS  IN  NAVFILES  IF  STORAGE  OPTION  IS  ON. 
C- 

IF(MIN(5).EQ.2HPS.0R.MIN(5).EQ.2HSP.0R.  MI N( 5)  .EQ  .  1HS ) CALL  BQ     0068 
*(ISYD,IHMS,SEMIMA,OECCEN,ORBINC ,MEANOM , PERHEL , ASNODE )  0069 

IOUT=l  0070 

C 

C      OUTPUT  RESULTANT  ORBIT  PARAMETERS  TO  THE  TERMINAL  CRT. 
C 

CALL  ORBOUT(IOUT,ISYD,IHMS,SEMIMA,OECCEN,ORBINC,M£a!\0;i,?L'RHEL,     0071 
*     ASNODE)  0072 

I0UT=2  0073 

C 

C      IF  PRINTER  OPTION  IS  ON,  ALSO  OUTPUT  ORBIT  PARAMETERS  TO  PRINTER. 
C 

IF(MIN(5).E0.2HPS.OR.MIN(5)  .  EQ.  2HSP  .OR  .  MI  N(  5)  .EQ . IflP  )CALL  O'RBOUT   0074 

* (IOUT, IS YD, IHMS, SEMI  MA, OECCEN .ORBIN C ,MEANOM , PERHEL , ASNODE)         0075 

RETURN  0076 

END  0077 

SIZE     2221     04255 
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SUBROUTINE  OEBOUT  ( IOUT  ,  IS  YD ,  IHMS  .SEMI.'IA  ,OECCEN  ,ORBINC  ,MEANOM, 

*   PERHEL,ASNODE) 

THE  SUBROUTINE  ORBOUT  TRANSMITS  THE  NEWLY  COMPUTED  ORBIT 
PARAMETERS  TO  THE  TERMINAL  CRT  OR  TO  THE  PRINTER  DEPENDING  ON  THE 
VALUE  OF  IOUT,   IOUT=l  MEANS  TERMINAL.   IOUT  =2  MEANS  PRINTER. 

INPUTS:   ALL  PARAMETERS 

OUTPUTS:  NO  NEW  PARAMETERS  OR  VALUES  ARE  RETURNED  TO  CALLING 

PROGRAM.   THE  ORBITAL  PARAMETERS  ARE  TRANSMITED  TO  THE 

OUTPUT  DEVICE  SPECIFIED  BY  IOUT. 

REAL  MEANOM 

INTEGER  M0UT(132) ,IARRAY(8) 

CALL  YDDMY(ISYD,ID,IM,IY) 

IARRAY(l)=(IY-1900)*10000+lM*100+ID 

IARRAY(2)=IHMS 

IARRAY(3)=IROUND(SEMIMA*100.0) 

IARRAY(4)=IROUND(OECCEN*1000000.0) 

IARRAY(5)=IROUND(ORBINC*1000.0) 

IARRAY(6)-IRQUND(MEANOM*1000.0) 

I ARRAY ( 7 )=IROUND( PERHEL*1000 .0 ) 

IAPRAY(8)=IROUND(ASNODE*1000.0) 

TP  TRANSMITS  A  HOLLERETH  FIELD  TO  THE  OUTPUT  DEVICE,  CRT  OR 
PRINTER 


0078 
0079 


0060 
0081 
0082 
0063 
0084 
0085 
0086 
0087 
0088 
0089 
0090 


CALL  TP(I0UT,132H   UPGORB  OUTPUT  RESULTS 


CALL  TP(I0UT,132H  KEPLERIAN  ORBIT  PARAMETERS 


MVCHAR  MOVEW  A  HOLLERETH  FIELD  INTO  THE  ARRAY  MOUT 

FIRST  BLANK  OUT  FIELD 
CALL  MVCHAR ( 'BLA' , 'BLA' , MOUT ,1,132) 
WRITE  FIELD  LABELS 
CALL  MVCHAR( 
*'    ETIMY    ETIMH   SEMIMA    ECCEN   ORBINC    MEANA  PERIGEE 
*,1,  MOUT, 1,64) 
CALL  TP(IOUT,MOUT) 
BLANK  FIELD  AGAIN 

CALL  MVCHAR ( 'BLA' , 'BLA ', MOUT , 1 , 132 ) 
CALL  TP(IOUT,MOUT) 
DO  10  1=1,8 

ENCODE  ORBIT  PARAMETERS  IN  INTEGER  HOLLERETH  FIELD 
CALL  MVCHAR(IARRAY(I ) , 'INT ', MOUT , 1*8, 1 ) 
10  CONTINUE 

TRANSMIT  HOLLERETH  FIELD  TO  OUTPUT  DEVICE 

CALL  TP(IOUT,MOUT) 

RETURN 


0091 
0092 
)0093 
0094 
0095 
)  0096 


0097 

0098 

ASNODE'0099 

0100 

0101 

0102 
0103 
0104 

0105 
0106 

0107 
0108 
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END  0109 

SIZE     400     00620 
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SUBROUTINE  DQ( ISYD ,IHMS , SEMI  MA ,OECCEN , ORB  I NC , MEANOM .PERHEL ,ASNODE ) 0110 

C 

C     DO  TRANSMITS  THE  ORBITAL  PARAMETERS  TO  THE  NAVIGATIONAL  FILES 

C 

C      INPUTS:   ORBITAL  PARAMETERS 

C     OUTPUTS:   THE  ORBIT  PARAMETERS  ARE  TRANSMITED  TO  THE  NAVIGATIONAL 

C  FILES  BY  NAVFIL. 

C 

REAL  MEANOM  0111 

INTEGER  MESONE(10)  ,MESTWO(10)  0112 

DATA  MESONE/'NAVFIL', 4,7*0/  0113 

DATA  MESTWO/'NAVFIL',5,7*0/  0114 

CALL  YDDMY(ISYD,ID,IM,IY)  0115 

I ETYMD=(IY-1900)* 10000+1 M*100+ID  0116 

ISEMI=IROUND(SEMIMA*100.0)  0117 

I OECC=IROUND(OECC EN* 1000000.0)  0118 

IASND=IROUND(ASNODE*1000.0)  0119 

I  01 NC  =  I ROUND (ORBING* 1000.0)  0120 

IPERH=IROUND(PERHEL*1000.0)  0121 

I  MEAN- IROUND( MEANOM* 1000.0)  0122 

MES0NE(4)=ISYD  0123 

MES0NE(6)=1  0124 

MES0NE(7)=IETYMD  0125 

MES0NE(8)=IHMS  0126 

MES0NE(9)=ISEMI  0127 

MES0NE(5)=MES0NE(9)  0128 

MESONE(10)=IOECC  0129 

CALL  SQ(MESONE)  0130 

MESTW0(4)=ISYD  0131 

MESTW0(6)=1  0132 

MESTW0(7)=I0INC  0133 

MESTW0(8)=IMEAN  0134 

MESTW0(9)=IPERH  0135 

MESTW0(5)=MESTW0(9)  0136 

MESTWO(10)=IASND  0137 

CALL  SQ(MESTWO)  0138 

RETURN  0139 

END  0140 

SIZE      139    00213 
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SUBROUTINE  GTLMBT ( ISYD ,NUMDAY ,HHHH , NUMLND , PTI MLN , XLINLM, XELELM f 
*     LCODE,ILNDBT,BETA,T) 

C 

C     GTLMBT  GETS  BETA  COUNTS  WITH  CALLS  TO  GETBET  AND  CHECKS  TO  SEE 

C      IF  THE  BETA  COUNT  PAIRS  OCCUR  AT  A  PICTURE  TIME  OF  ANY  OF  THE 

C      LANDMARKS.   IF  SO,  A  BETA  VALUE  FOR  THE  LANDMARK  IS  FOUND  BY 

C      INTERPOLATION  OR  EXTRAPOLATION.   ALSO,  THE  LANDMARK  CODE  IS 

C     EXAMINED  TO  DISCARD  STARS. 

C      INPUTS:   ISYD.   SATELLITE  ID.  YEAR  DAY 

C  NUMDAY.   THE  NUMBER  OF  DAYS  OF  LANDMARK  MEASUREMENTS  TO 

C  CHECK. 

C  HHHH.   THE  BEGINNING  HOUR  AND  ENDING  HOUR  TO  CONSIDER 

C  LANDMARK  MEASUREMENTS. 

C  NUMLND:   THE  NUMBER  OF  LANDMARKS  THAT  CAN  BE  HANDLED. 

C     OUTPUTS:   PTIMLM.  PICTURE  START  TIME  OF  LANDMARK 

C  XLINLM.   LINE  NUMBER  OF  LANDMARK 

C  XELELM.  ELEMENT  NUMBER  OF  LANDMARK 

C  LCODE.   LANDMARK  CODE 

C  ILNDBT.   LANDMARK-BETA  ASSOCIATION  FLAG.   EQUALS  1  WITH 

C  BETA.   EQUALS  -1  WITHOUT  BETA. 

C  T.      LANDMARK    TIME. 

C 

IMPLICIT  DOUBLE  PRECISION  (Q) 

DOUBLE  PRECISION  BETA( 1 ) ,RDPBT , BETDIF, T( 1 ) ,BETCMP 
REAL  PTIMLN(l) ,XLI NLM( 1 ), XELELM ( 1 ) 

INTEGER  HHHH  0147 

INTEGER  LCODE(l) ,ILNDBT(1)  0148 

INTEGER  NOMTIM(100)  0149 

INTEGER  ISCAN1 ( 100 ) , ISTIM1 ( 100  )  , ISMIL1 ( 100)  ,IBET1 (100)  0150 

INTEGER  ISCAN2(100),ISTIM2(100) ,ISMIL2(100) ,IBET2(100)  0151 

COMMON/NAVCOM/NAVDAY,LlNTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET 
1IMH,SEMIMA,0ECCEN,0RBINC,PERHEL,ASN0DE,N0PCLN,DECLIN,RASCEN,PICLIN 
2, PRERAT,PREDIR, PITCH, YAW, ROLL, SKEW 
COMMON/NAVINI/  0152 

1  EMEGA,AB,ASQ,BSQ,R,RSQ,  0153 

2  RDPDG,  0154 

3  NUMSEN,TOTLIN,RADLIN,  0155 

4  TOTELE,RADELE,PICELE,  0156 

5  CPITCH,CYAW,CROLL,  0157 

6  PSKEW,  0158 

7  RFACT,ROASIN,TMPSCL,  0159 

8  B11,B12,B13,B21,B22,B23,B31,B32,B33,  0160 

9  GAMMA, GAMDOT,  0161 
A  R0TM11,R0TM13,R0TM21,R0TM23,R0TM31,R0TM33,  0162 
B   PICTIM,XREF  0163 

DATA  MAXNUM/100/,NEGBET/3l44960/,RDPBT/9.989292881D-7/  0164 

DATA  BETCMP/3144960.0D0/  0165 
DATA  TWPI/6. 28318530/ 

ISAT=ISYD/100000  0166 

IYR=MOD(ISYD, 100000) /1000  0167 

IDAY=MOD(ISYD,1000)  0168 
C      COMPUTE  BEGINNING  AND  END  TIME  FOR  TIME  INTERVAL  CHECK 

BTIME=HHHH/100  0169 
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ETIME=BTIME+MOD(HHHH,100)  0170 

DO  10  I=1,NUMLND  0171 

10  ILNDBT(I)=-1  0172 
C     THE  BETAS  ARE  FETCHED  ONE  DAY  AT  A  TIME. 

DO  50  I=1,NUMDAY  0173 

JSYD=ISAT*100000+IYR*1000+IDAY  0174 

CALL  GETBET(JSYD,MAXNUM,NUM,NOMTIM,  0175 

*  ISCAN1,ISTIM1,ISMIL1,IBET1,  0176 

*  ISCAN2,ISTIM2,ISMIL2,IBET2)  0177 
IF(NUM.EQ.0)  GO  TO  40  0178 
DO  30  J=1,NUM  0179 

C     COMPUTE  BETA  TIME  USING  THE  I  INDEX  TO  COUNT  THE  HOURS  OF  EACH  DAY 

PTIME=FTIME(NOMTIM(J))+(I-1)*24.0  0180 

C      IF  TIME  LIES  OUTSIDE  SET  INTERVAL  FOR  ORBIT  COMPUTATION,  EXIT. 

IF(PTIME.LT.BTIME.OR.PTIME.GT.ETIME)GO  TO  30  0181 

C      THE  BETA  PAIRS  ARE  TAKEN  ONE  AT  A  TI^E  USING  THE  INDEX  J  AND  ALL 

C     THE  LANDMARKS  ARE  LOOKED  THROUGH  USING  THE  INDEX  K. 

C 

DO  25  K=1,NUMLND  0182 

C      CHECK  EACH  LANDMARK  OR  STAR  FOR  MATCH. 
C     IF  IT  IS  A  STAR,  DISCARD. 

IF(MOD(LCODE(K),100) .GE.5)G0  TO  25  0183 

C      IF  THE  TIMES  DON'T  MATCH,  DISCARD. 

IF(PTIME.NE.PTIMLN(K))GO  TO  25  0164 

C      IF  LANDMARK  HAD  PRIOR  MATCH,  DISCARD. 

IF(ILNDBT(K).EQ.l)GO  TO  25  0185 

C     OTHERWISE  LINEARLY  INTERPOLATE  OR  EXTRAPOLATE  MATCHING  BETA 
C     AND  TIME  VALUES  AND  SET  MATCH  ARRAY  ILNDBT  EQUAL  TO  ONE  AT  INDEX  K 

IF(IBET1(J).LT.0)  IBET1(J)=(IBET1(J).AND. '37777777) +NEGBET        0186 

IF(IBET2(J).LT.0)  IBET2 (J )=( IBET2 (J ) .AND . '37777777 ) +NEGBET        018? 
BETDIF=IBET2(J)-IBET1(J)  0188 

IF(BETDIF.GT.BETCMP)BETDIF=BETDIF-2.0D0*BETCMP  0189 

IF(BETDIF.LT.-BETCMP)BETDIF=BETDIF+2.0D0*BETCMP  0190 

ISCNLN=(IR0UND(XLINLM(K))+3)/8-ISCANl(J)+l 

QTIME1=FTIME(ISTIM1(J) ) +( 1-1 )*24. 0D0 

QTIME2=FTIME(ISTIM2(J) )+( 1-1 )*24 .0D0 

QDTIME=QTIME2-QTIME1 

QTM1ML=ISMIL1(J)/360000.0D0 

QTM2ML=ISMIL2(J)/360000.0D0 

QDTIME=(QTM2ML-QTM1ML)+QDTIME 

QTIME1=QTIME1+QTM1ML 

QTMPSL=QDTIME/(ISCAN2(J)-ISCAN1(J)) 

QBTPSL=BETDIF/(ISCAN2(J)-ISCAN1(J)) 

QPCTLN=(XELELM(K)-1.0)/(TOTELE-1.0)*(DEGELE/360.0) 

T(K)=QTIME1+QTMPSL*(ISCNLN+QPCTLN) 

BETA(K)=(IBET1(J)+QBTPSL*(ISCNLN))*RDPBT 

T(K)=T(K)+QTMPSL*BETA(K)/TWPI 
C     LANDMARK  POINTER  ARRAY  WHICH  IS  USED  FOR  RESIDUAL  COMPUTATION. 

ILNDBT (K)=l  0197 

25  CONTINUE  0198 

30  CONTINUE  0199 

C      INCREMENT  DAY  WHILE  ACCOUNTING  FOR  LEAPYEAR 
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IDAY=IDAY+1 

IF(IDAY.LE.365)G0  TO  40 

IF(IDAY.EQ.366.AND.MOD(IYR,4).EQ.0)GO  TO  40 

IDAY=1 

IYR=IYR+1 

IYR=MOD(IYR,100) 
40  CONTINUE 
50  CONTINUE 

RETURN 

END 

SIZE     1317     02445 
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SUBROUTINE  ANGORB( ISYD .NUMLND .MEANOM ,PTIME , XLAT ,XLON , 
*     XLIN.XELE.ILNDBT.BETA.V.T.NV.LNDPNT.IORBPL) 


0205 


WRITTEN  04/17/78  BY  DENNIS  PHILLIPS  TO  SET  UP  GEOMETRY  TO  FIND 

ORBITAL  POSITION  ANGLE,  I.  E.  TRUE  ANAMOLY+OFFSET  POSITON, 

IN  THE  ORBITAL  PLANE  OF  THE  SATELLITE. 

MODIFICATION  BY  DENNIS  PHILLIPS  ON  MAY  2,  1979  TO  ACCOUNT  FOR  THE  0209 

VARIATION  OF  THE  SATELLITE'S  HEIGHT.  0210 

MODIFICATION  07/17/78  BY  DENNIS  PHILLIPS  TO  DETERMINE 

ORIENTATION  OF  THE  ORBIT  PLANE  FROM  LANDMARK  MEASUREMENTS  IN      0212 

ADDITION  TO  THE  ALONG-PLANE  ORBIT  PARAMETER  DETERMINATION.        0213 

INPUTS:   ISYD.  SATELLITE  ID.   YEAR  DAY 

NUMLND.   NUMBER  OF  LANDMARKS 

MEANOM   MEAN  ANAMOLY  OF  SATELLITE  POSITION  AT   EPIC 

PTIME.   PICTURE  START  TIME  OF  LANDMARK 

XLAT.   LATITUDE  OF  LANDMARK 

XLON.   LONGITUDE  OF  LANDMARK 

XLIN.   LINE  NUMBER  OF  LANDMARK 

XELE.   ELEMENT  NUMBER  OF  LANDMARK 

ILNDBT.   LANDMARK-BETA  CORRESPONDENCE  FLAG 

BETA.   BETA  ANGULAR  SWEEP  TO  PICTURE  START  IN  RADIANS 

T.   TIME  OF  LANDMARK  MEASUREMENT 

IORBPL.   ORBIT  PLANE  CALCULATION  FLAG 
OUTPUTS:   V.   THE  TRUE  ANAMOLIES+OFFSET 

NV.   NUMBER  OF  LANDMARKS  WITH  CORRESPONDING  BETA  VALUES 
LNDPNT.   INDEX  POINTING  TO  LANDMARKS  WITH  MATCHING  BETAS 
ORBINC.   ORBIT  INCLINATION  IN  NAVCOM  COMMON 
ASNODE.   ASCENDING  NODE  IN  NAVCOM  COMMON. 

INTEGER  LNDPNT(l)  0214 

INTEGER  ILNDBT(l)  0215 
REAL  MEANOM 

REAL  PTIME(l),XLAT(l),XLON(l),XLlN(l),XELE(l)  0216 

DIMENSION  ESTST1(100),ESTST2(100),ESTST3(100)  0217 

DOUBLE  PRECISION  XCOR.YCOR  0218 

DOUBLE  PRECISION  BETA( 1 ) ,V (1  )  ,T (l )  0219 

DOUBLE  PRECISION  COSB.SINB.BETAR.TEMP , YCSNST .XCSNST  .XNORM  0220 

DOUBLE  PRECISION  SAMTIM  0221 
DOUBLE  PRECISION  UX.UY.UZ.PS I 

DOUBLE  PRECISION  XSN.YSN.ZSN  0223 

COMMON/MINCOM/UX(80) ,UY(80) ,UZ(80) ,PSI (80) ,NP  0225 
COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET0226 
lIMH,SEMIMA,OECCENtORBINC,PERHEL,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN0227 

2 .PRERAT .PREDIR , PITCH ,  YAW .ROLL .SKEW  0228 

COMMON/NAVINI/  0229 

1  EMEGA,AB,ASQ,BSQ,R,RSQ,  0230 

2  RDPDG  0231 

3  NUMSEN.TOTLIN.RADLIN,  0232 

4  TOTELE.RADELE.PICELE,  0233 

5  CPITCH.CYAW.CROLL,  0234 

6  PSKEW,  0235 

7  RFACT.ROASIN.TMPSCL,  0236 
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8  D11,D12,D13,D21,D22,D23,D31,D32,D33,  0237 

9  GAMMA ,GAMDOT,  0238 
A  R0TMlltR0TM13,R0TM21,R0TM23tR0TM31,R0TM33,  0239 
B  PICTIM.XREF  0240 

COMMON/BETCOM/IAJUSTrIBTCON,NEGBET,ISEANG  0241 

DATA  EH/42165.0/  0242 

DATA  PI/3.14159265/  0243 

DATA  IENT/1/ 

C     SET  UP  PARAMETERS  BEFORE  LANDMARK  LOOP 

IYD=MOD(ISYD, 100000)  0244 

NV=0  0245 

C      INIALIZE  ORBIT  PARAMETERS'  IE  THE  SEMI  MAJOR  AXIS  IS  AVAILABLE. 
NAVDAY=NAVDAT-1 

IF(SEMIMA.GT.42000.0)CALL  STVEC (0 .0 ,MEANOM ,XC ,YC , ZC ) 
NAVDAY=NAVDAY+1 
SENANG=FLALO(ISEANG)*RDPDG  0247 

C      THIS  IS  THE  LANDMARK  LOOP.   FOR  EACH  LANDMARK  WITH  MATCHING  BETA 

C     VALUES  AN  ANGULAR  POSITION  OF  THE  SATELLITE  IN  ITS  ORBITAL  PLANE 

C      IS  COMPUTED. 

DO  2  I=1,NUMLND  0263 

IF(ILNDBT(I)  .EQ.-DGO  TO  2  0264 

C      COUNT  MATCHES  BETWEEN  LANDMARKS  AND  BETAS 

NV=NV+1  0265 

C     TIME  VALUES  CORRESPONDING  TO  ORBITAL  ANGULAR  POSITIONS  ARE  ONLY 

C      SORTED  ON  THE  FIRST  ENTRY.   OTHERWISE  THESE  VALUES  GET  SHUFFLED 

C     AND  LOST. 

IF(IENT.EQ.1)T(NV)=T(I) 

TIME=T(NV) 

SAMTIM=T(NV) 

C      LANDMARK  POINTER  ARRAY  WHICH  IS  USED  FOR  RESIDULA  COMPUTATION. 

LNDPNT(NV)=I  0269 

C      PRECESS  ATTITUDE 

CALL  ATPREC(TIME,DEC,RAS) 

C     SET  UP  SATELLITE  CENTERED  COORDINATE  SYSTEM  DETERMINED  BY 

C     DECLINATION  AND  RIGHT  ASCENSION 
DEC=DEC*RDPDG 
RAS=RAS*RDPDG 

SR=SIN(RAS)  0250 

CR=COS(RAS)  0251 

SD=SIN(DEC)  0252 

CD=COS(DEC)  0253 

X1=-SR  0254 

Y1=CR  0255 

Z1=0.0  0256 

X2=-SD*CR  0257 

Y2=-SD*SR  0258 

Z2=CD  0259 

X3=CD*CR  0260 

Y3=CD*SR  0261 

Z3=SD  0262 

C     FIND  VECTOR  FROM  EARTH'S  CENTER  TO  LANDMARK  IN  INERTIAL  COODINATES 

C     TO  LANDMARK. 

JTIME=ITIME(TIME)  0271 
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RA-RAERAC(IYD,JTIME,XL0N(I))*RDPDG 

YLAT=XLAT(I)*RDPDG 

YLAT=GE0LAT(YLAT,1) 

TANLAT=TAN(YLAT)**2 

RR=SQRT((l.0+TANLAT)/(BSO+ASQ*TANLAT) )*AB 

X=COS(RA)*COS(YLAT)*RR 

Y=SIN(RA)*COS(YLAT)*RR 

Z=SIN(YLAT)*RR 
C     GET  VECTOR  FROM  EARTH  TO  SUN 

CALL  SUNVEC ( IYD ,SAMT IM , XS N ,YSN  ,ZSN ) 
C      SET  SATELLITE  VECTOR  TO  ZERO  FOR  FIRST  ITERATION 

XSAT=0.0 

YSAT=0.0 

ZSAT=0.0 

H=HH 

IF(SEMIMA.LT.1.0)GO  TO  1 
C      OTHERWISE  COMPUTE  THE  SATELLITE  VECTOR  IN  INERTIAL  COORDINATES 
C     AND  USE  HEIGHT  TO  IMPROVE  POSITION  ESTIMATE. 

CALL  STVEC(TIME,MEANOM,XSAT,YSAT,ZSAT) 

H  =  SORT (XSAT**2+YSAT**2+ZSAT**2  ) 
1  CONTINUE 
C      COMPUTE  X,  Y  COORDINATES  OF  THE  VECTOR  FROM  THE  SATELLITE  TO  THE 
C     SUN  IN  OUR  SATELLITE  CENTERED  COORDINATE  SYSTEM. 

XCSNST=X1*(-XSN-XSAT)+Y1*(-YSN-YSAT) 

YCSNST=X2*(-XSN-XSAT)+Y2*(-YSN-YSAT)+Z2*(-ZSN-ZSAT) 

THE  NEXT  SIX  STATEMENTS  YIELD  THE  POINTING  DIRECTION  OF  THE 

SPIN  SCAN  CAMERA  IN  A  COORDINATE  SYSTEM  WITH  THE  Z-AXIS  COINCIDIN( 

WITH  THE  VECTOR  OPPOSITE  THE  SPIN  AXIS  AND  THE  X-AXIS 

PERPENDICULAR  TO  THE  Z-AXIS  AT  AN  ANGLE  EQUAL  TO  SENANG  FROM  THE 

SUN  PULSE  DETECTOR. 

YLIN=(PICLIN-XLIN(I))*RADLIN 

CLIN=COS(YLIN) 

SLIN=SIN(YLIN) 

U=R0TM11*CLIN+R0TM13*SLIN 

V¥=R0TM21*CLIN+R0TM23*SLIN 

W=R0TM31*CLIN+R0TM33#SLIN 

THE  ANGULAR  SWEEP  DISTANCE  OF  THE  SPIN  SCAN  CAMERA  FROM  THE  SUN 

AT  THE  TIME  WHEN  THIS  PARTICULAR  LANDMARK  IS  BEING  VIEWED.   THE 

(VV.U)  TERM  ACCOUNTS  FOR  THE  ROLL  AND  YAW  MISALIGNMENT  EFFECTS. 

BETAR=BETA(I)-SENANG+XELE(I)*RADELE-ATAN2(VV,U) 

THE  NEXT  8  STATEMENTS  COMPUTE  A  UNIT  VECTOR  POINTING  AT  THE 

LANDMARK  IN  THE  SATELLITE  CENTERED  COORDINATE  SYSTEM. 

COSB=DCOS(BETAR) 

SINB=DSIN(BETAR) 

TEMP=COSB*XCSNST+SINB*YCSNST 

YCSNST=-SINB*XCSNST+COSB*YCSNST 

XCSNST=TEMP 

XNORM=DSQRT( (1 .0D0-W**2) /(XCSNST**2+YCSNST**2) ) 

XCSNST=XCSNST*XNORM 

YCSNST=YCSNST*XNORM 
C     FINALLY,  WE  HAVE  THE  POINTING  DIRECTION  OF  THE  SPIN  SCAN  CAMERA 
C      IN  INERTIAL  COORDINATES. 

XLNDST=-XCSNST*X1-YCSNST*X2-W*X3 


0272 
0273 
0274 
0275 
0276 
0277 
0278 
0279 

0280 

0281 
0282 
0283 

0284 
0285 


0287 
0292 


0293 
0294 


0295 
0296 
0297 
0298 
0299 
0300 


0301 


0302 
0303 
0304 
0305 
0306 
0307 
0308 
0309 


0310 


18.050575  PAGE       4 

YLNDST=-XCSNST*Y1-YCSNST*Y2-W*Y3  0311 

ZLNDST=-XCSNST*Z1-YCSN*T*Z2-W*Z3  0312 

C      COMPUTE  SLANT  RANGE  (w-xitNG)  TO  LANDMARK  FROM  SATELLITE. 

COSANG=XLNDST*X+YLNDST*Y+ZINDST*Z  0313 

SLTRNG=-C0SANG+SQRT(H**2-RR**2+C0SANG**2)  0314 

C  THE    INTERSECTION   OF   THE   RAY    EXTENDING    FROM   THE    LANDMARK   WITH      A 

C  SPHERE    CENTERED   AT   THE    EARTH'S    CENTER    AND   RADIUS    EQUAL   TO    H. 

ESTST1(NV)=X+XLNDST*SLTRNG  0315 

ESTST2(NV)=Y+YLNDST*SLTRNG  0316 

ESTST3(NV )=Z+ZLNDST*SLTRNG  0317 

IF(IORBPL.NE.l)GO  TO  2  0318 

C      THIS  VECTOR  IS  NORMALIZED  AND  THE  ANGLES  PSI(I)  ARE  SET  WHEN  THE 
C      ORIENTATION  OF  THE  ORBITAL  PLANE  IS  TO  BE  CALCULATED. 

YNORM=1.0/SQRT(ESTST1(NV)**2+ESTST2(NV)**2+ESTST3(NV)**2)         0319 
UX(NV)=ESTSTl(NV)*YNORM  0320 

UY(NV)=ESTST2(NV)*YN0RM  0321 

UZ(NV)=ESTST3(NV)*YN0RM  0322 

PSI(NV)=PI/2.0 

2  CONTINUE  0324 
IENT=0 

C      BRANCH  IF  THE  INCLINATION  AND  ASCENDING  NODE  ARE  NOT  TO  BE 
C      CALCULATED. 

IF(IORBPL.NE.l)GO  TO  3  0325 

NP=NV  0326 

NOPCLN=0 
C      CHANGE  ON  MARCH  22,  1980  BY  DENNIS  PHILLIPS  TO  REPLACE  CALL  TO 
C      MINMIZ,  THE  OLD  WAY  TO  COMPUTE  ATTITUDE.   INSTEAD  FINDAT,  A 
C      SUBROUTINE  USING  A  VARIATION  OF  NEWTON'S  METHOD,  IS  CALLED  TO 
C      COMPUTE  THE  ORIENTATION  OF  THE  ORBITAL  PLANE. 

CALL  FINDAT(NOPCLN  ,RADLIN ,ORBNC ,ASNDE , PCLN ) 

ORBINC=90.0+ORBNC 

ASNODE=ASNDE+270.0 

IF(ASNODE.GE.360.0)ASNODE=ASNODE-360.0 

3  CONTINUE  0338 
ASN=ASNODE*RDPDG  0339 
CASN=COS(ASN)  0340 
SASN=SIN(ASN)  0341 
XINC=ORBINC*RDPDG  0342 
CINC=COS(XINC)  0343 
SINC=SIN(XINC)  0344 
U1=CASN  0345 
Vl-SASN  0346 
W 1=0.0  034? 
U2=-SASN*CINC  0348 
V2=CASN*CINC  0349 
W2=SINC                                                         0350 

C 

C      COMPUTE  THE  ORBITAL  ANGLE  POSITIONS  OF  THE  SATELLITE  AROUND 

C      THE  ORBIT  PERPENDICULAR.   CARE  IS  TAKEN  TO  AVOID  A  WRONG  MULTIPLE 

C     OF  2*PI. 

C 

DO  10  1=1, NV  0351 

XC0R=U1*ESTST1 ( I ) +V1*ESTST2 ( I )+Wl*ESTST3(  I  )  0352 
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10 


YC0R=U2*ESTST1(I)+V2*ESTST2(I )+W2*ESTST3( I ) 

V(I)=DATAN2(YC0R,XC0R) 

IF(I.NE.1)G0  TO  5 

ANG2^V(l) 

TIME=T(1) 

GO  TO  7 

CONTINUE 

K=IROUND((T(I)-T(l))/24.0+(ANG2-V(I ))/(2.0*PI)) 

V(I)=V(I)+2.0*PI*K 

CONTINUE 

CONTINUE 

RETURN 

END 


0353 
0354 
0355 
0356 
0357 
0358 
0359 
0360 
0361 
0362 
0363 
0364 
0365 


SIZE 


1456 
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SUBROUTINE  STVEC ( SAMTIM .MEANOM ,X , Y , Z ) 
C   COPIED  0460  BY  DENNIS  PHILLIPS  FROM  SATVEC. 

C      STVEC  COMPUTES  THE  INERTIAL  POSITION  (X,  Y,  Z)  AT  TIME  SAMTIM. 
C      INPUTS:   SAMTIM   ZULU  TIME  IN  HOURS  FOR  WHICH  POSITION  IS  REQUIRED 
C  MEANOM   MEAN  ANAMOLY  OF  SATELLITE  POSITION  AT   EPIC 

C      OUTPUTS:   X,  Y,  Z    THE  INERTIAL  POSITION  OF  THE  SATELLITE  AT  ZULU 
C      TIME  SAMTIM. 

DOUBLE  PRECISION  RDPDG ,RE,GRACON .DIFTI M ,EACAN1 .ECANOM , XMANOM 

DOUBLE  PRECISION  DABS ,DSQRT ,DS IN ,DCOS 

REAL  MEANOM 

COMMON /NAVCOM/NAVDAY, LI NTOT.DEGLIN, IELTOT, DEGELE, SP INRA, IETIMY, I ET 
1IMH,SEMIMA,0ECCEN,0RBINC,PERIGE,ASN0DE,N0PCLN,DECLIN,RASCEN,PICLIN 
2,  PRERAT,PREDIR, PITCH, YAW, ROLL, SKEW 

DATA  NAVSAV/0/ 

DATA  GRACON,RE/.07436574D0,6378.388/ 

DATA  RDPDG/. 0174532925D0/ 
C      STVEC  WILL  RE-INITIALIZE  ORBIT  PARAMETERS  IF  NAVDAY  IS  CHANGED. 

IF(NAVDAY.EQ.NAVSAV)  GO  TO  10 

0=RDPDG*ORBINC 

P=RDPDG*PERIGE 

A=RDPDG*ASNODE 

SO-SIN(O) 

CO=COS(0) 

SP=SIN(P)*SEMIMA 

CP=COS(P)*SEMIMA 

SA=SIN(A) 

CA=COS(A) 

PX=CP*CA-SP*SA*CO 

PY=CP*SA+SP*CA*CO 

PZ=SP*SO 

QX=-SP*CA-CP*SA*CO 

QY=-SP*SA+CP*CA*CO 

QZ=CP*SO 

SROME2=SQRT(1.0-OECCEN)*SQRT(1.0+OECCEN) 

XMMC=GRAC0N*RE*DSQRT(RE/SEMIMA)/SEMIMA 

EPSILN=1.0E-8 
C    COMPUTE  THE  DIFFERENCE  IN  TIMES  OF  ( NAVDAY , 0 )-( IETIMY , IETIMH) . 
C      THIS  DIFFERENCE  IS  COMPUTED  IN  MINUTES. 

IEY=MOD( IETIMY /10 00, 100) 

INY=MOD(NAVDAY/1000,100) 

IYRMN=0 

IF(IEY.EQ.INY)GO  TO  5 

IF(INY.LT.IEY)GO  TO  3 

IE=INY-1 

DO  2  I=IEY,IE 

IYRMN=IYRMN+365*1440 

IF(MOD(I,4).EQ.0)IYRMN=IYRMN+1440 

2  CONTINUE 
GO  TO  5 

3  CONTINUE 
IE=IEY-1 

DO  4  I=INY,IE 
IYRMN=IYRMN-365*1440 
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IF(MOD(I,4).EQ.0)IYRMN=IYRMN-1440 

4  CONTINUE 

5  CONTINUE 
IED=MOD(IETIMY,1000) 
IND=MOD(NAVDAY,1000) 
IYRMN=IYRMN+(IND-IED)*1440 
TE=FLOAT(IYRMN)-FTIME(IETIMH)*60.0 

10  DIFTIM=SAMTIM*60.0+TE 

XMANOM=XMMC*DIFTIM+MEANOM*RDPDG 

ECANMl=XMANOM 

DO  20  1=1,20 

ECANOM=XMANOM+OECCEN*DSIN(ECANMl) 

IF(DABS(ECANQM-ECANM1) .LT .EPSILN )GO  TO  30 
20  ECANMl=ECANOM 
30  XQMEGA=DCQS(ECANOM)-OECCEN 

Y0MEGA=SR0ME2*DSIN(ECAN0M) 

X=XOMEGA*PX+YOMEGA*QX 

Y=XOMEGA*PY+YOMEGA*QY 

Z=XOMEGA*PZ+YOMEGA*QZ 

RETURN 

END 

SIZE     345    00531 
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SUBROUTINE  ORBPR( V ,T ,N ,MEANOM , ISEMI  ,EPTLME) 
C      WRITEEN  BY  DENNIS  PHILLIPS  OF  SCIENTIFIC  PROGRAMMING  AND  APPLIED 
C     MATHEMATICS   INC. 

C      ORBPR  CONVERTS  ANGULAR  MEASUREMENTS  (TRUE  ANAMOLY+OFFSET)  IN  THE  0 
C      ORBPR  CONVERTS  ANGULAR  MEASUREMENTS  (TRUE  ANAMOLY+OFFSET)  IN  THE 
C     ORBITAL  PLANE  WHICH  HAVE  ASSOCIATED  TIME  TAGS  INTO  KEPLERIAN 
C     ORBITAL  ELEMENTS. 

C      INPUTS:   V.  TRUE  ANAMOLY+OFFSET  ANGLES. 
C  T.   TIME  TAGS  OF  ANGLES 

C  N.   NUMBER  OF  ANGLES. 

C  ISEMI.   NUMBER  OF  KEPLERIAN  ORBITAL  PARAMETERS  TO  3E 

C  COMPUTED. 

C  EPTIME.   EPIC  TIME  FOR  KEPLERIAN  ORBITAL  ELEMENTS. 

C      OUTPUTS:   MEANOM.   MEAN  ANAMOLY  OF  ORBIT  AT  EPIC  TIME. 
C  SEMIMA.   ORBIT'S  SEMIMAJOR  AXIS  IN  NAVCOM  COMMON  BLOCK. 

C  OECCEN.   ORBIT'S  ECCENTRICITY  IN  NAVCOM  COMMON  BLOCK 

C  PERIGEE.   ORBIT'S  PERIGEE  IN  NAVCOM  COMMON  BLOCK. 

DOUBLE  PRECISION  T(1),V(1) 

REAL  MEANOM 

COMMON/NAVCOM/NAVDAY.LINTOT ,DEGLIN , IELTOT .DEGELE , SPINRA , IETIMY , I ET 
1IMH,SEMIMA,OECCEN,ORBINC,PERIGE,ASNODE,NOPCLN,DECLIN,RASCEN,PICLIN 
2,PRERAT,PREDIR, PITCH, YAW, ROLL, SKEW 

DATA  GRACON ,RE/ .07436574D0 ,6378.388/ 

DATA  PI, RDPDG/3. 14159265,. 01745329252/ 
C      COMPLETELY  REWRITTEN  SEPTEMBER  25,  1979  BY  DENNIS  PHILLIPS  OF  SPAA 
C      M  TO  CLARIFY   THE  CODING  METHODS  AND  TO  CONFORM  TO  A  STRUCTURED 
C      PROGRAMMING  PHILOSOPHY.    ALSO,  THE  METHOD  FOR  COMPUTINGE  THE 
C      SEMIMAJOR  AXIS  IS  MUCH  IMPROVED. 

IF(ISEMI.NE.4)GO  TO  20 
C      COMPUTE  ALL  FOUR  ALONG-TRACK  PARAMETERS  BY  LINEAR  REGRESSION. 

CALL  F0URCF(C1,C2,C3,C4,V,T,N) 

SEMIMA=(GRACON*C2*60.0)**(2.0/3.0)*PE 

CALL  ORBTHR ( CI, C2,C3,C4, EPTIME, MEANOM) 

RETURN 
20  IF(ISEMI.NE.3)G0  TO  30 
C     COMPUTE  ORBIT'S  MEAN  ANAMOLY,  ECCENTRICITY,  AND  ARGUMENT  OR  PERIGE 

C2=(SEMIMA/RE)**(3.0/2.0)/(60.0*GRACON) 

CALL  THRC0F(C1,C2,C3,C4,V,T,N) 

CALL  ORBTHR ( CI, C2,C3,C4, EPTIME, MEANOM) 

RETURN 
30  IF(ISEMI.NE.2)G0  TO  40 
C      COMPUTE  ORBIT'S  SEMIMAJOR  AXIS  AND  MEAN  ANAMOLY  WHILE  ADJUSTING 
C      FOR  THE  PERTURBATIVE  EFFECTS  OF  A  NONZERO  ECCENTRICITY. 

C3=0.0 

C4=0.0 

DO  36  J=l,4 

SME=0.0 

SMESQ=0.0 

SMT=0.0 

SMET=0.0 

XN=FLOAT(N) 

DO  35  1=1, N 

ANG=V(I) 
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TT=T( I )-C3*SIN ( ANG)-C4*C0S (ANG) 
SME=SME+ANG 
SMESQ=SMESQ+ANG**2 
SMT=SMT+TT 

35  SMET=SMET+ANG*TT 

DET=1 .0/(XN*SMESQ-SME*SME ) 
C2=(SMET*XN-SME*SMT)*DET 
C3=-2.0*OECCEN*C2*COS(PERIGE*RDPDG) 
C4=2.0*C2*OECCEN*SIN(PERIGE*RDPDG) 

36  CONTINUE 
C1=(SMT*SMESQ-SMET*SME)*DET 
SEMIMA=(GRACON*C2*60.0)**(2.0/3.0)*RE 
MEANOM= (EPTIME-C1 ) /( C2*RDPDG )-PERIGE 
MEANOM=AMOD(MEANOM, 360.0) 
IF(MEANOM.LT.0.0)MEANOM=MEANOM+360.0 
RETURN 

40  IF(ISEMI.NE.1)RETURN 

COMPUTE  ORBIT'S  MEAN  ANAMOLY. 

C2=(SEMIMA/RE)**(3.0/2.0)/(60.0*GRACON) 

C3=-2.0*OECCEN*C2*COS(PERIGE*RDPDG) 

C4=2.0*C2*OECCEN*SIN(PERIGE*RDPDG) 

CALL  0NEPAR(C2,C3,C4,V ,T , N ,EPT 1MB, PERIGE.MEANOM) 

RETURN 

END 

SIZE     321     00501 
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SUBROUTINE  FOURCF ( CI ,C2  ,C3 ,C4  ,V ,T ,N ) 
C      WRITTEN  BY  DENNIS  PHILLIPS  OF  SCIENTIFIC  PROGRAMMING  AND  APPLIED 
C     MATHEMATICS,  INC. 

C      FOURCF  USES  LINEAR  REGRESSION  TO  MINIMIZE  THE  SUM  OVER  I  OF 
C      (CI  +  C2*V(I)  +  C3*C0S(V(I))  +  C4*SIN(V(I))  -  T(I))**2. 
C      INPUTS:   TRUE  ANAMOLY+OFFSET  ANGULAR  POSITIONS. 
C  T.   ANGULAR  TIME  TAGS. 

C  N.   NUMBER  OF  ANGULAR  POSITIONS. 

C      OUTPUTS:   CI,  C2,  C3  AND  C4   LINEAR  PLUS  SINE  WAVE  CONSTANTS  OF 
C  ORBITAL  MOTION. 

DOUBLE  PRECISION  A (4 ,4 ) ,B (4 ,4 ) ,C (4 ) ,D (4 ) , V ( 1 )  ,T (l  ) 

DO  5  1=1,4 

JS  =  I 

D(I)=0.0 

DO  5  J=JS,4 
5  A(I,J)=0.0 
C      SUM  REGRESSION  TERMS  INTO  THE  APPROPIATE  ELEMENTS  OF  THE  MATRIX  A 

XK>N 

A(1,1)=XK 

DO    10    1=1, N 

ANG=V(I) 

TT=T(I) 

CANG=COS 

SANG=SIN 

A(1,2)=A 

A(2,2)=A 

A(1,3)=A 

A(1,4)=A 

A(2,3)=A 

A(2,4)=A 

A(3,3)=A 

A(3,4)=A 

A(4,4)=A 

D(1)=D(1 

D(2)=D(2 

D(3)=D(3 

D(4)=D(4 
10  CONTINUE 
C      INVERT  A 

CALL  INVE(A,B) 

DO  20  1=1,4 

C(I)=0.0 

DO  20  K=l,4 
20  C(I)=C(I)+B(I,K)*D(K) 

C1=C(1) 

C2=C(2) 

C3=C(3) 

C4=C(4) 

RETURN 

END 


(ANG) 

(ANG) 

(1,2)+ANG 

(2,2)+ANG**2 

(1,3)+SANG 

(1,4)+CANG 

(2,3)+ANG*SANG 

(2,4)+ANG*CANG 

(3,3)+SANG**2 

(3,4)+SANG*CANG 

(4,4)+CANG**2 

)  +  TT 

)+TT*ANG 

)+SANG*TT 

)+TT*CANG 


SIZE 


266 


00412 
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SUBEOUTINE  INVE(A,B) 
C     WRITTEN  BY  DENNIS  PHILLIPS 
C     MATHEMATICS,  INC. 
C      INVERTS  FOUR  BY  FOUR  REGRES 
C     INPUTS:  FOUR  BY  FOUR  MATRIX 
C     OUTPUTS:   B,  THE  INVERSE  OF 
C     A  IS  PARTITIONED  INTO  FOUR 
C     RESULTANT  TWO  BY  TWO  MATRIX 
C     BY  USING  GAUSSIAN  ELIMINATI 
DOUBLE  PRECISION  A(4,4),B(4 
DOUBLE  PRECISION  DET ,TEMP 
DET=1.0D0/(A(1,1)*A(2,2)-A( 


B(l,l 
B(l,2 
B(2,l 
B(2t2 
C(l,l 
C(2,l 
C(l,2 
C(2,2 
B(3,l 
B(4,l 
B(3,2 
B(4,2 
D(l,l 
D(l,2 
D(2,l 
D(2,2 
DET  =  1 
B(3,3 
B(3,4 
B(4,3 
B(4,4 


TEMP=B(3,1) 


B(3,l 
B(4,l 


TEMP=B(3f2) 


B(3,2 

B(4,2 

B(l,3 

B(l,4 

B(2,3 

B(2,4 

B(l,l 

B(l,2 

B(2,l 

B(2,2 

RETURN 

END 


=A(2,2)*DET 

=-A(l,2)*DET 

=B(l,2) 

=A(1,1)*DET 

=B(1,1)*A(1,3)+B(1,2) 

=B(2,1)*A(1,3)+B(2,2) 

=B(1,1)*A(1,4)+B(1,2) 

=B(2,l)*A(l,4)+B(2,2) 

=-A(l,3)*B(l,l)-A(2,3 
=-A(l,4)*B(l,l)-A(2,4 

=-A(l,3)*B(l,2)-A(2,3 
=-A(l,4)*B(l,2)-A(2,4 
=A(3,3)-A(1,3)*C(1,1) 


OF  SCIENTIFIC  PROGRAMMING  AND  APPLIED 

SION    MATRIX. 

A. 

A. 
TWO    BY   TWO    SUBMATRICES    AND    THEN    THE 

CONSISTING    OF   SUBMATRICES    IS    INVERTED 
ON. 
,4),C(2,2),D(2,2) 

1.2)*A(1,2)) 


=A(3,4 
=A(3,4 
=A(4,4 
0D0/(D 
=D(2,2 


=-D(l,2)*DET 
=B(3,4! 


=D(1,1 


=B(3,3 
=B(4,3 


=B(3,3 
=B(4,3 
=B(3,1 
=B(4,1 
=B(3,2 
=B(4t2 
=B(1,1 
=B(1,2 
=B(1,2 
=B(2,2 


~A(l,3)*C(l,2) 
-A(1,4)*C(1,1) 
-A(1,4)*C(1,2) 
1,1)*D(2,2)-D( 
*DET 


*A(2,3) 

*A(2,3) 

*A(2,4) 

*A(2,4) 

)*B(2,1) 

)*B(2,1) 

)*B(2-,2) 

)*B(2,2) 

-A(2,3)*C(2,1) 

-A(2,3)*C(2,2) 

-A(2,4)*C(2,1) 

-A(2,4)*C(2,2) 

1,2)*D(2,1)) 


*DET 

*B(3,1)+B(3,4) 
*TEMP+B(4,4)*B 

*B(3,2)+B(3,4) 
*TEMP+B(4,4)*B 


*B(4,1) 

(4,1) 

*B(4,2) 
(4,2) 


-C(1,1)*B(3,1) 
-C(1,1)*B(3,2) 


-C(l,2)*B(4,i; 
-C(l,2)*B(4t2; 


-C(2,1)*B(3,2)-C(2,2)*B(4,2) 


SIZE 


271 


00417 


Fl 
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SUBROUTINE  ORBTHR(Cl ,C2 ,C3 ,C4 ,EPTIME, MEANOM) 

C      PULLED  OUT  OF  ORBPR  BY  DENNIS  PHILLIPS  ON  MAY  2,  1979  0607 

C      CHANGED  BY  DENNIS  PHILLIPS  ON  MAY  11,  1979  0608 

C     ANAMOLY  AND  EQUALS  ONE  WHEN  THE  COEFFICIENTS  ARE  COMPUTED  FROM  THE0610 

C     ECCENTRIC  ANAMOLY.  0611 

C      COMPUTE  ECCENTRICITY,  ARGUMENT  OF  PERIGEE  AND  MEAN  ANAMOLY.  0612 

C      COMPUTES  ECCENTRICITY,  ARGUMENT  OF  PERIGEE  AND  MEAN  ANAMOLY  FROM 
C      CI,  C2,  C3  AND  C4  COEFFICIENTS. 
C      INPUTS:   CI,  C2,  C3  AND  C4  COEFFICIENTS. 
C              EPTIME.   THE  EPIC  TIME  OF  ORBIT  COMPUTATION. 
C      OUTPUTS:  OECCEN .   ORBIT  ECCENTRICITY. 
C              PERIGE.'  ORBIT'S  PEPIGEE 
C              MEANOM.   ORBIT'S  MEAN  ANAMOLY. 
C 

REAL  MEANOM  0613 

COMMON /NAVCOM/NAVDAY ,LI NTOT .DEGLIN , IELTOT ,DEGELE , SPINRA , IETIMY , I ET0614 
1IMH,SEM IMA, OECCEN, ORBINC , PERIGE , ASNODE , NOPCLN .DECLIN .RASCEN ,PI CLIN0615 

2,PRERAT,PREDIR,PITCH,YAW,RCLL,SKEW  0616 

COMMON/NAVINI/  0617 

1  EMEGA,AB,ASC,BS0,R,RSO,  0616 

2  RDPDG,  0619 

3  NUMSEN,TOTLIN,RADLIN,  0620 

4  TOTELE,RADELE,PICELE,  0621 

5  CPITCH,CYAW,CROLL,  0622 

6  PSKEW,  0623 

7  RFACT,ROASIN,TMPSCL,  0624 

8  B11,B12,B13,B21,B22,B23,B31  ,B32,B33,  0625 

9  GAMMA, GAMDOT,  0626 
A  R0TM11,R0TM13,R0TM21 ,R0TM23 ,R0TM31 .R0TM33 ,  0627 
B   PICTIM,XREF  0628 

OECCEN=SGRT(C3**2+C4**2)/(2.0*C2^ 
PERIGE=0.0 

PERIGE=ATAN2(C4,-C3) /RDPDG 
PERIGE=AMOD( PERIGE ,360.0) 
IF(PERIGE.LT.0.0)PERIGE=PERIGE+360.0 
40  CONTINUE  0635 

MEANOM=(EPTIME-Cl )/( C2*RDPDG ) -PERIGE 

MEANOM=AMOD(MEANOM, 360.0)  0637 

IF (MEANOM. LT.0.0)MEANOM=MEANOM+360.0  0638 

RETURN  0639 

END  0640 

SIZE      78    00116 
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SUBROUTINE  THRCOF( CI ,C2,C3 ,C4,E ,T ,N )  0566 

WRITTEN  MAY  2,1980  BY  DENNIS  PHILLIPS  FOR  NESS. 

THRCOF  FITS  A  SINE  WAVE  PLUS  CONSTANT  TO  THE  CURVE  T (I  )-C2*E(I ) .  0568 
EFFECTIVELY  A  LEAST  SQUARES  FIT  IS  PERFORMED  TO  FIND  Cl9  C3  AND  C40569 

IN  THE  EXPRESSION  C1+C2*E( I ) +C3*SIN (E( I ) ) +C4*COS(E( I ) )=T ( I ) .  0570 
INPUTS:   E.   THE  TRUE  ANAMOLY  +  OFFSET  SATELLITE  POSITION  ARRAY. 

T.   THE  TIME  OF  THE  ESTIMATED  SATELLITE  POSITIONS. 

N.   THE  NUMBER  OF  ESTIMATED  SATELLITE  POSITIONS. 

C2.   THE  COEFFICIENT  DETERMINING  THE  SATELLITE'S  MEAN 
OUTPUTS:   CI,  C3  AND  C4.   THE  CONSTANT  PLUS  SINE  WAVE  COEFFICIENTS 

DOUBLE  PRECISION  E(1),T(1)  0571 

SMSN=0.0  0572 

SMCS=0.0  0573 

SMSNSQ=0.0  0574 

SMSNCS=0.0  0575 

SMCSSQ=0.0  0576 

SMCN=0.0  0577 

SMSNCN=0.0  0578 

SMCSCN=0.0  0579 
COLLECT  REGRESSION  COEFFICIENTS 

DO  10  1=1, N  0580 

ANG=E(I)  0581 

CN=T(I)-C2*ANG  0582 

SN=SIN(ANG)  0583 

CS=COS(ANG)  0584 

SMSN=SMSN+SN  0585 

SMCS=SMCS+CS  0586 

SMSNSQ=SMSNSQ+SN*SN  0587 

SMSNCS=SMSNCS+SN*CS  0588 

SMCSSQ=SMCSSQ+CS*CS  0589 

SMCN=SMCN+CN  0590 

SMSNCN=SMSNCN+SN*CN  0591 

SMCSCN=SMCSCN+CS*CN  0592 

10  CONTINUE  0593 

XN=N  0594 
INVERT  THRREE  BY  THREE  REGRESSION  MATRIX  USING  CRAMER'S  RULE. 

DET=XN*(SMSNSQ*SMCSSQ-SMSNCS**2)-SMSN*(SMSN*SMCSSQ-SMCS*SMSNCS)  0595 

*+SMCS*(SMSN*SMSNCS-SMCS*SMSNSQ)  0596 

DET=1.0/DET  0597 
C1=(SMCN*(SMSNSQ*SMCSSQ-SMSNCS**2)-SMSNCN*(SMSN*SMQSSQ-SMCS*SMSNCS0598 

*)+SMCSCN*(SMSN*SMSNCS-SMCS*SMSNSQ))*DET  0599 
C3=(XN*(SMSNCN*SMCSSQ-SMSNCS*SMCSCN )-SMSN* (SMCN*SMCSSQ-SMCS*SMCSCN06  00 

*)+SMCS*(SMCN*SMSNCS-SMCS*SMSNCN) )*DET  0601 
C4=(XN*(SMSNSQ*SMCSCN-SMSNCN*SMSNCS)-SMSN*(SMSN*SMCSCN-SMCN*SMSNCS0602 

*)+SMCS*(SMSN*SMSNCN-SMCN*SMSNSQ) )*DET  0603 

RETURN  0604 

END  0605 


SIZE 


250 


00372 


18.050575  PAGE  1 

SUBROUTINE  ONEPAR ( C2 ,C3 ,C4 fV ,T , N tEPTIME .PERIGE ,MEANOM ) 
C      WRITTEN  BY  DENNIS  PHILLIPS  OF  SCIENTIFIC  PROGRAMMING  AND  APPLIED 
C      MATHEMATICS,  INC. 

C      COMPUTE  THE  ORBIT'S  MEAN  ANAMOLY  BY  REGRESSION. 

C      INPUTS:   C2,  C3  AND  C4   LINEAR  AND  SINE  WAVE  CONSTANTS  OF  ORBITAL 
C  MOTION. 

C  V.   TRUE  ANAMOLY  +  OFFSET  ANGULAR  POSITIONS. 

C  T.   ANGULAR  TIME  TAGS. 

C  N.   NUMBER  OF  ANGULAR  POSITIONS. 

C  EPTIME.   EPIC  TIME  OF  DESIRED  MEAN  ANAMOLY  POSITION. 

C  PERIGE.   ARGUMENT  OF  PERIGEE. 

C     OUTPUTS:   MEANOM:  MEAN  ANAMOLY. 

DOUBLE  PRECISION  V(1),T(1) 

REAL  MEANOM 

DATA  RDPDG/. 01745329252/ 

SUM=0.0 

XN-N 

DO  10  1=1, N 

ANG=V(I) 
10  SUM=SUM+T(I)-C2*ANG-C3*SIN(ANG)-C4*C0S(ANG) 

C1=SUM/XN 

MEANOM=(EPTIME-Cl )/(C2*RDPDG )-PERIGE 

MEANOM=AMOD (MEANOM ,360 .0 ) 

IF ( MEANOM. LT.0.0)MEANOM=MEANOM+360.0 

RETURN 

END 

SIZE      100     00144 
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SUBROUTINE  ORBRES ( ISYD , IOUT , NUMLND, T,XLI N ,XELE , BETA , MEANOM 9 

*  LCODE,XLAT,XLON,LNDPNT)  036? 
C 

C      ADAPTED  10/05/78  BY  DENNIS  PHILLIPS  FROM  SUBROUTINE  RESIDU        0368 

C     OUTPUTS  LINE  AND  ELEMENT  RESIUDES  FROM  UPGORB  0369 

C 

C      INPUTS:   ISYD.   SATELLITE  ID.   YEAR   DAY 

C  IOUT.   OUTPUT  DEVICE.   ONE  MEANS  CRT.   TWO  MEANS  PRINTER. 

C  NUMLND.   THE  NUMBER  OF  LANDMARKS  TO  COMPUTE .RESIDUALS  FOR 

C  T.   THE  TIME  OF  EACH  LANDMARK  MEASUREMENT. 

C  XLIN.   THE  LINE  NUMBER  OF  EACH  LANDMARK  MEASUREMENT. 

C  XELE.   THE  ELEMENT  NUMBER  OF  EACH  LANDMARK  MEASUEEMNT. 

C  BETA.   THE  SWEEP  ANGLE  UP  TO  THE  BEGINNING  OF  THE  SCAN 

C  LINE  ON  WHICH  THE  LANDMARK  WAS  MEASURED. 

C  MEANOM   MEAN  ANAMOLY  OF  SATELLITE  POSITION  AT   EPIC 

C  LCODE.   THE  LANDMARK  CODE. 

C  XLAT.   THE  LATITUDE  OF  THE  MEASURED  LANDMARK. 

C  XLON.   THE  LONGITUDE  OF  THE  MEASURED  LANDMARK. 

C  LNDPNT.   AN  INDEX  ARRAY  POINTING  AT  LANDMARKS  WITH 

C  ASSOCIATED  BETA  COUNTS. 

C 

INTEGER  LCODE(l),LNDPNT(l)  0370 

REAL  MEANOM 

REAL  XLAT(l)fXLON(l),XLIN(l>,XELE(l)  0371 

DOUBLE  PRECISION  T(1),BETA(1)  0372 

INTEGER  LINE(44)  0373 

C0MM0N/NAVC0M/NAVDAY,LINT0T,DEGLIN,IELT0T,DEGELE,SPINRA,IETIMY9IET 
1IMH,SEMIMA,0ECCEN,0RBINC,PERIGE,ASN0DE,N0PCLN,DECLIN,RASCEN9PICLIN 
2, PRERAT,PREDIR, PITCH, YAW, ROLL, SKEW 
COMMON/NAVINI/ 

1  EMEGA,AB,ASQ,BSQ.R,RSQ, 

2  RDPDG, 

3  NUMSEN.TOTLIN.RADLIN, 

4  TOTELE,RADELE,PICELE, 

5  CPITCH,CYAW,CROLL, 

6  PSKEW, 

7  RFACT,ROASIN,TMPSCL, 

8  B11,B12,B13,B21,B22,B23,B31,B32,B33, 

9  GAMMA, GAMDOT, 

A   R0TM11,R0TM13,R0TM21,R0TM23,R0TM31,R0TM33, 
B   PICTIM,XREF 

IF(NUMLND.EQ.0)  GO  TO  91  3374 

C     RE-INITIALIZE  STVEC  WITH  NEW  ORBIT  PARAMETERS 

NAVDAY-NAVDAY-1 

CALL  STVEC(0.0, MEANOM, X,Y,Z) 

NAVDAY=NAVDAY+1 

CALL  STVEC(0.0, MEANOM, X,Y,Z) 
C 

C      TP  SENDS  A  HOLLERETH  FIELD  TO  THE  OUTPUT  DEVICE  IOUT. 
C 

CALL  TP(I0UT,132H  LANDMARK  RESIDUALS  FROM  UPGORB  0375 

*  0376 

*  0377 
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C      SEND  SECOND  HOLLERETH  FIELD 

CALL  TP(I0UT,132H   N  SSYYDDD   HHMMSS    LCODE   LINDIF   ELEDIF  LANDL0376 
*AT   LANDLON  0379 

*  )0380 

INAV=1  0362 

IYD=MOD(ISYD, 100000)  0383 

DO  10  I=1,NUMLND  0384 

C 

C      MVCHAR  MOVES  A  HOLLERETH  FIELD  INTO  THE  ARRAY  LINE. 
C 
C      FIRST  BLANK  OUT  THE  FIELD 

CALL  MVCHARCBLA',  'BLA',  LINE,  1,132)  0385 

C     PLACE  INTEGER  HOLLERETH  FIELD  IN  ARRAY 

CALL  MVCHAR(I, 'INT', LINE, 3,1)  0386 

C     PLACE  INTEGER  HOLLERETH  FIELD  IN  ARRAY 

CALL  MVCHAR(ISYD, 'INT', LINE, 11,1)  0387 

TIME=T(I)  0388 

JTIME-ITIME(TIME)  0389 

C     PLACE  INTEGER  HOLLERETH  FIELD  IN  ARRAY 

CALL  MVCHAR (JTIME, 'INT', LINE, 19,1)  0390 

K=LNDPNT(I)  0391 

C     PLACE  INTEGER  HOLLERETH  FIELD  IN  ARRAY 

CALL  MVCHAR(LCODE(K) , 'INT', LINE, 27,1)  0392 

ILINE=(XLIN(K)-1.0)/8.0  +  1.0 

PTIME  =  T( I )-TMPSCL*FLOAT( I  LINE) 
C      TRANSFORM  LANDMARK'S  LATITUDE  AND  LONGITUDE  TO  IMAGE  COORDINATES 
C      USING  THE  NEW  ORBIT  PARAMETERS. 

CALL  CALRES(IYD,T(I) .BETA (K ) ,MEANOM ,XLAT(K) ,XLON(K) ,YLIN,YELE) 
C      CALCULATE  THE  DISTANCE  BETWEEN  THE  CALCULATED  AND  MEASURED 
C      IMAGE  POSITIONS. 

DLIN=XLIN(K)-YLIN  0396 

DELE=XELE(K)-YELE  0397 

C      TRANSFER  RESIDUALS  INTO  A  FLOATING  HOLLERETH  FIELD. 

CALL  FF(8,2,DLINtLINE,27)  0400 

CALL  FF(8f2, DELE, LINE, 35)  0401 

C      TRANSFER  THE  LATITUDE  AND  LONGITUDE  OF  LANDMARK  INTO  AN  INTEGER 

CALL  MVCHAR(ILALO(XLAT(K) ), 'INT', LINE, 51,1)  0402 

C      HOLLERETH  FILE. 

CALL  MVCHAR(ILALO(XLON(K) ),  'INT',  LINE, 60,1)  0403 

C      OUTPUT  GENERATED  HOLLERETH  LINE. 

CALL  TP(IOUT,LINE)  0404 

10  CONTINUE  0405 

RETURN  0406 

91  CALL  TQMESCNO  LANDMARKS  FOR  SSYYDDD:  $',  ISYD )  0407 

RETURN  0408 

END  0409 

SIZE     419     00643 
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SUBROUTINE  CALRES ( IYD.T ,BETA ,  MEANOM ,XLAT ,XLCN ,XLI N ,XELE) 
C     WRITTEN  BY  DENNIS  PHILLIPS  OF  SCIENTIFIC  PROGRAMMING  AND  APPLIED 
C     MATHEMATICS,  INC.   CALRES  COMPUTES  THE  LINE  AND  ELEMENT  POSITION 
C     FOR  A  GIVEN  (XLAT ,  XLON)  USING  THE  BETA  SWEEP  ANGLE. 
C     INPUTS:   IYD.   YEAR  DAY 

C  T.   TIME  OF  LANDMARK  SCAN  LINE. 

C  BETA.   SWEEP  ANGLE  ASSOCIATED  WITH  LANDMARK  SCAN  LINE. 

C  MEANOM   MEAN  ANAMOLY  OF  SATELLITE  POSITION  AT   EPIC 

C  XLAT.   THE  LATITUDE  OF  THE  EARTH  LANDMARK. 

C  XLON.   THE  LONGITUDE  OF  THE  EARTH  LANDMARK. 

C     OUTPUTS:   XLIN.   THE  COMPUTED  LINE  POSITION. 
C  XELE.   THE  COMPUTED  ELEMENT  POSITION. 

REAL  MEANOM 

DOUBLE  PRECISION  T ,XSN , YSN , ZSN ,XSNST , YSNST ,XST  , YST ,ANG  ,TWPI , BETA   0411 

COMMON/NAVCOM/NAVDAY,LINTOT,DEGLIN,IELTOT,DEGELE,SPINRA,IETIMY,IET0412 
1IMH,SEMIMA.,0ECCEN,0RBINC  , PERHEL  .ASNODE  ,  NOPCLN  .DEC  LIN  ,RASCEN  ,PICLI  N0413 
2, PRERAT.PREDIR, PITCH, YAW, ROLL, SKEW  0414 

COMMON/NAVINI/  0415 

1  EMEGA,AB,ASQ,BSQ,R,RSQ,  0416 

2  RDPDG,  0417 

3  NUMSEN,TOTLIN,RADLIN,  0418 

4  TOTELE,RADELE,PICELE,  0419 

5  CPITCH.CYAW.CROLL,  0420 

6  PSKEW,  0421 

7  RFACT,ROASIN,TMPSCL,  0422 

8  D11,D12,D13,D21,D22,D23,D31,D32,D33,  0423 

9  GAMA,GAMDOT,  0424 
A  R0TM11,R0TM13,R0TM21,R0TM23,R0TM31,R0TM33,  0425 
B   PICTIM,XREF                                                   0426 

COMMON/BETCOM/IAJUST,IBTCON,NEGBET,ISEANG  0427 

DATA  TWPI/6.28318531D0/  0428 

TIME=T  0430 

JTIME=ITIME(TIME)  0431 

C     GET  VECTOR  TO  THE  SUN  IN  INERTIAL  COORDINATES 

CALL  SUNVEC(IYD,T,XSN,YSN,ZSN)  0432 

C     GET  VECTOR  TO  THE  SATELLITE  IN  INERTIAL  COORDINATES 

CALL  STVEC( TIME, MEANOM, XS AT, YSAT.ZS AT) 
C     COMPUTE  THE  X,  Y,  Z  VECTOR  POINTING  TO  THE  LANDMARK  IN  INERTIAL 
C     COORDINATES  FROM  THE  EARTH'S  CENTER. 

JTIME-ITIME(TIME) 

YLON=RAERAC(IYD,JTIME,XLON)*RDPDG 

TDIF=TIME-FTIME(JTIME) 

YLON=YLON+TDIF*TWP 1/24.0 

YLAT=XLAT*RDPDG 

YLAT=GEOLAT(YLAT,l) 

TANLAT=TAN(YLAT)**2 

RR=SQRT( (1.0+TANLAT)/(BSQ+ASQ*TANLAT) )*AB 

XE=COS(YLON)*COS(YLAT)*RR 

YE=SIN(YLON)*COS(YLAT)*RR 

ZE=SIN(YLAT)*RR 

SENANG=FLALO(ISEANG)*RDPDG  0439 

C      FETCH  CURRENT  ATTITUDE  STATE  AND  SETUP  SATELLITE  CENTERED 
C      COORDINATE  SYSTEM  WITH  Z-AXIS  OPPOSITE  THE  SPIN  AXIS. 
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CALL  ATPREC(TIME,DEC ,RAS^ 

DEC=DEC*RDPDG 

RAS=RAS*RDPDG 

SR-SIN(RAS)  0442 

CR-COS(RAS)  0443 

SD=SIN(DEC)  0444 

CD=COS(DEC)  0445 

X1=-SR  0446 

Y1=CR  0447 

Z1=0.0  0448 

X2=-SD*CR  0449 

Y2=-SD*SR  0450 

Z2=CD  0451 

X3=CD*CR 

Y3=CD*SR 

Z3=SD 
C      FIND  LINE  NUMBER 

C0SANG=X3*(XE-XSAT)+Y3*(YE-YSAT)+Z3*(ZE-ZSAT) 

COSANG=COSANG/SQRT( (XE-XSAT )**2+ ( YE-YSAT)**2+ ( ZE-ZSAT )**2 ) 

THETA-ATAN2 (R0TM31 ,R0TM33 ) 

PHI=ASIN(C0SANG/SQRT(R0TM13**2+R0TM33**2) ) 

YLIN=PHI-THETA 

XLIN=PICLIN-YLIN/RADLIN 
C      FIND  ELEMENT  NUMBER 

C      FIND  COORDINATES  OF  THE  VECTOR  FROM  THE  SATELLITE  TO  THE  SUN 
C      IN  THE  SATELLITE  COORDINATE  SYSTEM 

XSNST=X1*(-XSN-XSAT)+Y1*(-YSN-YSAT)  0452 

YSNST=X2*(-XSN-XSAT)+Y2*(-YSN-YSAT)+Z2*(-ZSN-ZSAT)  0453 

C      FIND  THE  COORDINATES  OF  THE  VECTOR  FROM  THE  SATELLITE  TO  THE  EARTH 
C      LANDMARK  IN  THE  SATELLITE  COORDINATE  SYSTEM. 

XST=X1* (XE-XSAT )+Y 1* ( YE-YSAT ) 

YST=X2*(XE-XSAT)+Y2*(YE-YSAT)+Z2*(ZE-ZSAT) 
C      CORRECT  FOR  ROLL  AND  SKEW 

CLIN=COS(YLIN) 

SLIN=SIN(YLIN) 

U=R0TM11*CLIN+R0TM13*SLIN 

V=R0TM21*CLIN+R0TM23*SLIN 

ANG=DATAN2(YSNST,XSNST)-DATAN2(YST,XST)-BETA+SENANG+ATAN2(V,U) 
C      NORMALIZE  THIS  ANGLE  TO  LIE  BETWEEN  -PI/2.0  AND  PI/2.0 

ANG-DMOD(ANG,TWPI ) 

IFUNG.GT.  (TWPI/2.0D0)  )  ANG=ANG-TWPI 

IF(ANG.LT.  (-TWPI/2.0D0)  )    ANG=ANG+TW'PI 

XELE=ANG/RADELE 

RETURN  0465 

END  0466 


SIZE     436    00664 
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SUEROUTINF  EARLOC ( ISC , FRTIME, DATF , SLINF ,FLEM , ELAT , ELONG ) 
C 
C 

C        PURPOSE: 

C  FARLOC  COMPUTES  THE  GEODETIC  LATITUDE  (ELAT)  AND  LONGITUDE 

C  (FLOMG)  CORRESPONDING  TO  THE  GIVEN  LINF  (SLINF)  AND  ELEMENT 

C  (ELEM)  COORDINATES  OF  THE  VISSR  IMAGE  FRAME  STARTING  AT  FRTIME, 

C  EARLOC  PFRFORMS  THE  INVERSE  TRANSFORM  OF  IMGLOC. 

C 
C 

C        PROGRAMMER:  LARRY  HAM  BRICK,  NESS /OS  I 
C 

C        ORIGINATION  DATE:  FEBRUARY  18,1976 
C 
C 

C        ARGUMENTS: 
C 

C  THE  UNITS  OF  ELAT  AND  ELONG  ARE  DECIMAL  DEGREES  WITH  ELAT 

C  POSITIVE  NORTH  AND  FLONG  POSITIVE  EAST. 

C 

C  THE  UNITS  OF  SLINF  AND  ELEM  ARE  SCAN  LINES  AND  IR  SAMPLES 

C  RESPECTIVELY.  THEIR  FRACTIONAL  PARTS  SPECIFY  HIGHER 

C  RESOLUTIONS. 

C 

C  NOTE:  THE  RANGE  FOR  ELEM  IS  0.5  TO  3822.5.  THE  IR  SAMPLE 

C  NUMBER  IS  OBTAINED  BY  ROUNDING-OFF  ELEM.  THE  VISIBLE 

C  SAMPLE  NUMBER  IS  OBTAINED  BY  ROUNDING-OFF  ((ELEM-0.5) 

C  *4  +  0.5). 

C 

C  THE  RANGF  FOR  SLINE  IS  0.5  TO  1821.5.  THE  IR  (SCAN) 

C  LINE  NUMBER  IS  OBTAINED  BY  ROUNDING-OFF  SLINE.  THE 

C  VISIBLE  LINE  NUMBER  IS  OBTAINED  BY  ROUNDING-OFF 

C  ( (SLINE-0.5)*8  +  0.5). 

C 

C  ISC  DESIGNATES  THE  SATFLLITE:  1  FOR  EAST  GOES,  2  FOR  WEST. 

C 

C  FRTIME  IS  THE  GMT  IN  SECONDS  OF  THE  START  OF  THE  IMAGE  FRAME. 

C 

C  DATE  SPECIFIES  THE  YEAR  AND  JULIAN  DAY  AS  YYDDD. 

C 
C 
C 

C        REFERENCES: 
C 

C  THE  MATHEMATICAL  BASIS  FOR  THIS  ROUTINE  IS  THE  PAPER  BY 

C  DENNIS  PHILLIPS  AND  C . T .MOTTERSHEAD .  THE  REPORT  BY  WESTING- 

C  HOUSF  (CONTRACT  NAS 5-23582 .DATED  JULY, 1977)  IS  MORE  DIRECTLY 

C  LINKED  TO  THE  FORTRAN  CODE  AND  GIVES  DETAILED  DEFINITIONS 

C  AND  ILLUSTRATIONS  OF  THE  ORBIT/ATTITUDE  PARAMETERS  AVAILABLE 

C  IN  THE  VISSR  DOCUMENTATION  BLOCK. 

C 

C  EARLOC  IS  NOT  IN  THE  MOST  EFFICIENT  COMPUTATIONAL  FORM; 

C  EMPASIS  HERE  IS  ON  CLARITY  OF  THE  GEOMETRY. 


18.050575  PAGE       2 

C 
C 

DIMENSION  PHI(6),S(3,3) 

COMMON  /  GRNICH  /  GRA1 ,GRA2 

COMMON  /  TIM  /  TIMF1,DATE1,D 

COMMON  /  SATATT  /  SPRA1 , SPRA2, SPDC1 , SPDC2 

COMMON  /  SATPOS  /  CX ( 11 ) ,C Y ( 11  )  ,CZ ( 1 1) 

COMMON  /  ALNANG  /  ZETA,RHO,ETA 

COMMON  /  SPNRAT  /  SPPER 
C 

C  GRA1  GREENICH  ANGLE  IN  DEGREES  AT  TIME  TIME1 

C  GRA2  GREENWICH  ANGLE  IN  DEGREES  AT  TIMF  TIME1  +  D 

C  TIME1  GMT  IN  SECONDS  OF  EPOCH  FOR  NAVIGATION  PARAMETERS 

C  D  PERIOD  IN  SECONDS  OVER  WHICH  TIME  IS  NORMALIZED 

C  DATF1  DATE  AS  YYDDD,YEAR  AND  JULIAN  DAY, OF  EPOCH  FOR  NAV .  PAR 

C  SPRA-,SPDC-  RIGHT  ASCENSION  AND  DECLINATION  OF  POSITIVE  SPIN 

C  AXIS  IN  DEGREES  AT  TIMES:  (l)TIMEl  AND  (2)  TIME1 

C  PLUS  D 

C  CX(I ),CY(I) ,CZ(I ),I=1  TO  11   CHEBYCHEV  COEFFICIENTS  REPRESENT 

C  -ING  THE  SATELLITE  POSITION  IN  KM  BEGINNING  AT  EPOCH 

C  (TIMED  IN  EARTH  CFNTFRED  INERTIAL  COORDINATES 

C  NORMALIZED  WRT  TIME  OVER  D 

C  RHO  IS  ROLL  OR  ELEMFNT  BIAS  OF  THE  VISSR  IN  DFGRFFS 

C  ZETA  IS  PITCH  OR  LINE  BIAS  OF  THE  VISSR  IN  DEGREES 

C  ETA  IS  YAW  OR  SKEW  OF  THE  VISSR  IN  DEGREES 

C  SPPFR  SPIN  PERIOD  OF  SATFLLITF  IN  SECONDS 

C 

DATA  ASB,ACQANG, TOTSMP, SCENCA,SCENCS, TOTSCL  /  6378.144,6356.759, 
1    9.1875,3822.0,45.0,4096.0,1821.0  / 
C 

C  A:  EQATORIAL  RADIUS  OF  OBLATE  ELLIPSOID  EARTH   IN  KM 

C  B:  POLAR  RADIUS  OF  OBLATE  ILLIPSOID  FARTH  IN  KM 

C  ACQANG:  DATA  ACQUISTION  ANGLE  OF  VISSR  IN  DEGREES 

C  TOTSMP:  TOTAL  IR  SAMPLES  ACQUIRED  IN  SCAM  LINI 

C  SCENCA:  VISSR  SCAN  ENCODER  CHARACTERISTIC  ANGLE  IN  DEGREES 

Z  SCFNCS:  NUMBER  OF  VISSR  SCAN  ENCODER  POSITIONS 

C  TOTSCL:  TOTAL  NUMBER  OF  SCAN  LINES 

C 

PI=3. 1415926535 
C 

C      COMPUTE  SOME  CONSTANT  PARAMETERS 
C 

E=(A**2  -  B**2)/B**2 

BSAS=  B**2/A**2 

AMUF=  (2.0*ACQANG/TOTSMP )*PI/180.0 

AMUL=    (SCENCA/SCENCS)*PI/180.0 

CE=  (TOTSMP  +  1.0)/2.0 

CL=  (TOTSCL  +  1.0)/2.0 
C 

C      COMPUTE  VISSR  ALIGNMENT  ANGLES 
C 

PITCH    =    ZETA*PI/180.0 

YAW=   ETA*PI/180.0 
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ROLL  =  RHO*PI/180.0 
C 

C  CHECK  THAT  ELEM  AND  SLIN?  ARE  IN  LIMITS 

C 

IF(  (ELFM.GT.3822.5hOR.  ( ELEM .LT .0 . 5)  )    GO  TO  10 

IF( (SLINE.GT.1821.5).0R. (SLINE .LT .0 .5) )  GO  TO  10 

GO  TO  15 
10  PRINT  30 
30   F0RMAT(41H   LINE  OR  ELEMENT  NUMBER  IS  OUT  OF  LIMITS  ) 

ELAT=999.9 

FL0NG=999.9 

GO  TO  20 
15    CONTINUE 
C 

C  COMPUTE  POLAR  COORDINATES  OF  THE  UNIT  VIEW  VECTOR  IN 

C  THF  VISSR  COORDINATE  SYSTEM 

C 

AZM=AMUE*(ELEM-CE) 

ELV  «  AMUL*(CL-SLINE^ 
C 

C  TRANSFORM  THE  UNIT  VFCTOR  TO  SATELLITE  INFRTIAL  COORDINATES 

C 

VSl=COS (AZM+ROLL )*COS (PI TCH+ELV )-SI N (AZM+ROLL )*SI N ( YAW )*SI N (PI TCH 
1  +ELV) 

VS2=+SI M( AZM+ROLL ) *COS ( P ITCH+ELV ) +COS ( AZM+ROLL) -SIN ( YAW )*SIN( PITCH 
1  +ELV) 

VS3=-C0S(YAW)*SIN(PITCH+ELV} 
C 

C  DETERMINE  TIME  OF  THE  GIVFN  SAMPLE  (APPROX.  TO  THE  SCAN 

C  LINE  TIME) 

C 

T  =  FRTIME  +   AINT(SLINE  +0.5)*SPPER 
C 

C  TRANSFORM  THF  UNIT  VIEW  VECTOR  TO  EARTH  CENTERED  INERTIAL 

C  COORDINATES 

C 

CALL  SATCRD(T,S ) 

ITC1=  S(l,l)*VSl  +  S(2,1)*VS2  +  S(3,lv:sVS3 

VC2  =  S(1,2)*VS1  +  S(2,2)*VS2  +  S(3,2)*VS3 

VC3  =  S(1,3)*VS1  +  S(2,3)*VS2  +  S(3,3)*VS3 
C 

C  EXTEND  THE  UNIT  VECTOR  TO  INTERSECTION  tflTH  THE  EARTH 

C  ELLIPSOID 

C 

CALL  SATVFC(T,PCX,PCY,PCZ) 
F=(PCX*VC1+PCY*VC2+PCZ*VC3  +  F*PCZ*VC3 ) /( 1 .0+E*VC3**2 ) 

G=(PCX**2+PCY**2+PCZ**2-A**2  +  F*PCZ**2) /(l .0  +  E*VC3**2) 

Q  =  -F  -  SQRT(F**2-G) 
C 

C  COMPUTE  THE  VECTOR  FROM  THF  EARTH  CENTER  TO  THE  SAMPLE 

C  POINT  IN  EARTH  CENTERED  INFRTIAL  COORDINATES 

C 

RC1  =  PCX  +  Q*VC1 
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RC2  -  PCY  +  Q*VC2 

RC3  =  PCZ  +  Q*VC3 
C 

C  COMPUTE  THE  EARTH  GEODETIC  LATITUDE  AND  LONGITUDE  OE  SAMPLE 

C 

ELAT  =  ATAN( (1.0/BSAS  )*RC3/SQRT(RC1**2+RC2**2))  *180.0/PI 
C 

C  THE  LONGITUDE  IS  COMPUTED  UNDER  THF  CONVENTION  -180  TO  +180 

C  WITH  EAST  OF  GREENWICH  + 

C 

ANG  =  ATAN2(RC2,RC1  ) 
C 

C  COMPUTF  THE  GREENWICH  ANGLE 

C 

CALL  GRANG(T.W) 
C 

ELONG  =  ANG  -  W 

IF(ELONG.LE.-PI )  EL0NG=2 .0*PI+ELONG 

IF(ELONG.GT.PI )  FL0NG=EL0NG-2 . 0*PI 

ELONG  =  ELONG*180.0/PI 
20   CONTINUE 

RETURN 

END 

SIZE     499    00763 
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SUBROUTINE  IMGLOC ( ISC ,FRTI ME, DATE , ELAT , ELONG , SLINE , ELEM ) 
C 
C 

C        PURPOSE: 

C  IMGLOC  COMPUTES  THE  LINE  (SLINE)  AND  ELEMENT  (ELEM)  COOR- 

C  DINATFS.IN  THE  VISSR  IMAGE  FRAME  STARTING  AT  TIME  (FRTIME), 

C  CORRESPONDING  TO  THE  GIVEN  EARTH  GEODETIC  LATITUDE  (ELAT)  AND 

C  LONGITUDE  (ELONG). 

C  IMGLCC  PERFORMS  THE  INVERSE  TRANSFORM  OF  EARLOC. 

C 
C 

C         PROGRAMMER:  LARRY  HAMBRICK,  NESS/OSI 
C 

C        ORIGINATION  DATE:  FEBRUARY  18,1978 
C 
C 

C        ARGUMFNTS: 
C 

C  THT  UNITS  OF  FLAT  AND  ELONG  ARE  DECIMAL  DEGREES  WITH  ELAT 

C  POSITIVE  NORTH  AND  ELONG  POSITIVE  EAST. 

C 

C  THE  UNITS  OF  SLINE  AND  ELEM  ARE  SCAN  LINES  AND  IR  SAMPLES 

C  RESPECTIVELY.  THEIR  FRACTIONAL  PARTS  SPECIFY  HIGHER 

C  RESOLUTIONS. 

C 

C  NOTE:  THE  RANGE  FOR  ELEM  IS  0.5  TO  3822.5.  THE  IR  SAMPLE 

C  NUMBER  IS  OBTAINED  BY  ROUNDING-OFF  ELEM.  THE  VISIBLE 

C  SAMPLE  NUMBER  IS  OBTAINED  BY  ROUNDING-OFF  ((ELEM-0.5) 

C  *4  +  0.5). 

C 

C  THE  RANGE  FOR  SLINE  IS  0.5  TO  1821.5.  THE  IR  (SCAN) 

C  LINE  NUMBER  IS  OBTAINED  BY  ROUNDING--OFF  SLINE.  THE 

C  VISIBLE  LINE  NUMBER  IS  OBTAINED  BY  ROUNDING-OFF 

C  ((SLINE-0.5)*8  +  0.5) . 

C 

C  ISC  DESIGNATES  THE  SATELLITE:  1  FOR  EAST  GOES,  2  FOR  JEST. 

C 

C  FRTIME  IS  THE  GMT  IN  SECONDS  OF  THE  START  OF  THE  IMAGE  FRAME. 

C 

C  DATE  SPECIFIES  THE  YEAR  AND  JULIAN  DAY  AS  YYDDD. 

C 
C 

C        REFERENCES: 
C 

C  THE  MATHEMATICAL  BASIS  FOR  THIS  ROUTINE  IS  THE  PAPER  BY 

C  DENNIS  PHILLIPS  AND  C .T .MOTTERSHEAD .  THE  REPORT  BY  tfESTING- 

C  HOUSE  (CONTRACT  NAS 5-23582 , DATED  JULY,19?7)  IS  MORE  DIRECTLY 

C  LINKED  TO  THE  FORTRAN  CODE  AND  GIVES  DETAILED  DEFINITIONS 

C  AND  ILLUSTRATIONS  OF  THE  ORBIT/ATTITUDE  PARAMETERS  AVAILABLE 

C  IN  THE  VISSR  DOCUMENTATION  BLOCK. 

C 

C  IMGLOC  IS  NOT  IN  THE  MOST  EFFICIENT  COMPUTATIONAL  FORM; 

C  EMPASIS  HERE  IS  ON  CLARITY  OF  THE  GEOMETRY. 
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C 
C 

DIMENSION  PHI (6),T(6),S(3,3) 

COMMON  /  GRNICH  /  GRA1,GRA2 

COMMON  /  TIM  /  TIMEl.DATEl  ,D 

COMMON  /  SATATT  /  SPRA1 , SPRA2, SPDC1 , SPDC2 

COMMON  /  SATPOS  /  CX( 11 ) , C Y ( 11)  ,CZ ( 1 1  ) 

COMMON  /  ALNANG  /  ZETA  ,RHO  ,  ETA. 

COMMON  /  SPNRAT  /  SPPER 
C 

C  GRA1  GREENICH  ANGLE  IN  DEGREES  AT  TIME  TIME1 

C  GRA2  GREENWICH  ANGLE  IN  DEGREES  AT  TIME  TIME1  +  D 

C  TIME1  GMT  IN  SECONDS  OF  EPOCH  FOR  NAVIGATION  PARAMETERS 

C  D  PERIOD  IN  SECONDS  OVER  WHICH  TIME  IS  NORMALIZED 

C  DATF1  DATE  AS  YYDDD.YEAR  AND  JULIAN  DAY, OF  EPOCH  FOR  NAV  .  PAR 

C  SPRA-.SPDC-  RIGHT  ASCENSION  AND  DECLINATION  OF  POSITIVE  SPIN 

C  AXIS  IN  DEGREES  AT  TIMES:  (l)TIMEl  AND  (2)  TIME1 

C  PLUS  D 

C  CX(I  ),CY(I),CZ(I) ,1=1  TO  11   CHEBYCHEV  COEFFICIENTS  REPRESENT 

C  -ING  THE  SATELLITE  POSITION  IN  KM  BEGINNING  AT  EPOCH 

C  (TIMED  IN  EARTH  CENTERFD  INERTIAL  COORDINATLS 

C  NORMALIZED  WRT  TIME  OVER  D 

C  RHO  IS  ROLL  OR  ELEMENT  BIAS  OF  THE  VISSR  IN  DEGREFS 

C  ZETA  IS  PITCH  OR  LINE  BIAS  OF  THE  VISSR  IN  DEGREES 

C  ETA  IS  YAW  OR  SKEW  OF  THE  VISSR  IN  DEGREES 

C  SPPER  SPIN  PERIOD  OF  SATELLITE  IN  SECONDS 

C 

DATA  STDFC, A, B.AC3ANG, TOTSMP, SCENCA.SCENCS, TOTSCL  / 
1   6.611,6378.144,6356.759,9.1375,3822.0,45.0,4096.0,1821.0  / 
C 

C  STDEC:RATIO  OF  NOMINAL  SATELLITE  RANGE  AND  EARTH  RADIUS 

C  A:  I0ATORIAL  RADIUS  OF  OBLATF  ELLIPSOID  EARTH   IN  KM 

C  B:  POLAR  RADIUS  OF  OBLATE  ELLIPSOID  EARTH  IN  KM 

DATA  ACQUISTION  ANGLE  OF  VISSR  IN  DEGREES 

TOTAL  IR  SAMPLES  ACQUIRED  IN  SCAN  LINE 

VISSR  SCAN  ENCODER  CHARACTERISTIC  ANGLE  IN  DEGREES 

NUMBER  OF  VISSR  SCAN  ENCODER  POSITIONS 

TOTAL  NUMBER  OF  SCAN  LINES 


C  ACQAMG 

C  TOTSMP 

C  SCENCA 

C  SCENCS 

C  TOTSCL 


C 

PI-  3.1415926535 
C 

C      COMPUTE  SOME  CONSTANT  PARAMETERS 
C 

E=(A**2  -  B**2)/B**2 

BSAS=  B**2/A*#2 

AMUE=  (2.0*ACOANG/TCTSMP )*PI/180.0 

AMUL=  (SCENCA/SCENCS )*PI/180.0 

CE=  (TOTSMP  +  1.0V2.0 

CL=  (TOTSCL  +  1.0)/2.0 
C 

C       COMPUTE  VISSR  ALIGNMENT  ANGLES 
C 

PITCH  =  ZETA*PI/180.0 
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YAtf  =  ETA*PI/180.0 

ROLL  =  RHO*PI/180.0 
C 

C  CONVERT  GEODETIC  LATITUDE  AND  LONGITUDE  TO  RADIANS 

C 

PLAMDA=  ELAT*PI/180.0 

GLCNG=  ELONG*PI/180.0 
C 

C  COMPUTE  GEOCENTRIC  LATITUDE 

C 

GLAMDA=ATAN(BSAS*TAN(PLAMDA)) 
C 

C  ESTIMATE  THE  TIME  AT  WHICH  THE  POINT  WAS  SCANNED 

C 

PHI(1)=  ATAN(  1.0*SIN(PLAMDA)/(STDEC-COS(PLAMDA)))-PITCH 

T(1)=FRTIME-(AINT(PHI(1)/AMUL  +  0 .5 )-CL )*SPPER 
C 

C  THIS  LOOP  ITERATES  TO  REFINE  THE  TIME  ESTIMATE 

C 

DO  40  1=1,5 

TI=T(I ) 
C 

C  COMPUTE  THE  GREENWICH  ANGLE 

C 

CALL  GRANG(TI  ,W) 
C 

C  DETERMINE  THE  EARTH  POINT  VFCTOR  IN  EARTH  CENTERED  INERTIAL 

C  COORDINATES 

C 

XLl=COS(GLAMDA)*COS(GLONG+W) 

XL2=C0S(GLAMDA)*SIN(GL0NG+W) 

XL3=SIN(GLAMDA) 

R=A/SQRT(1.0  +  E*SIN(GLAMDA)**2) 

RC1=R*XL1 

RC2=R*XL2 

RC3=R*XL3 
C 

C  DETERMINE  THE  SATELLITE  VECTOR 

C 

CALL  SATVEC(TI  ,PX,PY,PZ) 
C 

C  CHECK  WHETHER  EARTH  POINT  IS  IN  VIEW 

C 

XLDP=XL1*PX  +  XL2*PY  +  XL3*PZ 

IF(XLDP.LT.A)  GO  TO  50 
C 

C  COMPUTE  VIEW  VFCTOR  IN  EARTH  CENTERED  INERTIAL  COORDINATES 

C 

VC1=RC1-PX 

VC2=RC2-PY 

VC3=RC3-PZ 
i 
C  TRANSFORM  THE  VIEW  VECTOR  TO  SATELLITE  INERTIAL  COORDINATES 
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C 

CALL  SATCRD(TI ,S) 

VS1=S(1,1)*VC1  +  S(l,2)*VC2  +  S(1,3)*VC3 
VS2=S(2,ll*VCl  +  S(2,2)*VC2  +  S(2,3)*VC3 
VS3=S(3,1)*VC1  +  S(3,2)*VC2  +  S(3,3)*VC3 
C 

C  DFTFRMINE  AZIMUTH  AND  ELFVATION  OF  EARTH  POINT  IN  SAT.  COORD. 

C 

SIGMA=ATAN(VS2/VS1) 
TANYAW=TAN(YAW) 
COSYAW=COS(YAW) 

XI=ATAN(VS3*TANYAW/SQRT(VS1**2+VS2**2-VS3**2*TANYa«/**2)) 
PHI (1+1 )=ATAN(-VS3/(C0SYAW*C0S(XI )*SQRT( VS1**2+VS2**2) ) )  -PITCH 
C 

C  CHECK  FOR  CONVERGENCE  OF  THE  TIME  ESTIMATE 

C 

DSLINE=(AINT(PKI(I+l)/AMUL  +  0  .  5  )-AI  NT  (  PHI  ( I  )/AMUL  +  0.5)) 
K  =  I+1 

IF(DSLINE.GE.1.0)  GO  TO  15 
GO  TO  25 
C 

C  CORRICT  THE  ESTIMATE  OF  TIME 

C 
15    T(I+1)=T(I)-DSLINE*SPPER 
40   CONTINUE 
C 

C  COMPUTE  THE  LINE  AND  ELEMENT 

C 

25   THETA=XI+SIGMA-ROLL 
ELEM=  CE  +  THETA/AMUE 
SLINE=CL-(PHI (K)/AMUL) 
C 

C  ROUGH  CHECKS  ON  REASONABLENESS  OF  RESULTS 

IF(FLEM.LT.0.5)  GO  TO  35 
IF(ELEM.LE.3822.5)  GO  TO  30 
35   PRINT  39 

39   FORMAT(17H   ELFMENT  ERROR   ) 
30    IF(SLINE.LT.0.5)  GO  TO  45 

IF(SLINE.LE.1821.5)  GO  TO  60 
45   PRINT  49 

49   F0RMAT(14H    LINE  ERROR   ) 
GO  TO  60 
50  PRINT  55 
55   FORMAT(18H   BEHIND  THE  EARTH) 
ELEM=0.0 
SLINF=0.0 
60   CONTINUE 
RETURN 
END 

SIZE      556    01066 
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SUBROUTINJ  SATCRD(T.S) 
C  COMPUTES  THF  BASF  VECTORS  (  THE  MATRIX  S)  OF  THE  SATELLITE 

C  COORDINATE  SYSTEM  IN  TERMS  OF  THF  EARTH  CENTERED  INFRTIAL 

C  COORDINATE  SYSTEM  AT  TIME  (T)  BASED  ON  THE  SATELLITE 

C  POSITION  AND  ATTITUDE 

COMMON  /  SATPOS  /  CX ( 11) ,C Y ( 1 1)  ,CZ ( 11 ) 

COMMON  /  SATATT  /  SPRA1 , SPRA2 , SPDC1 , SPDC2 

COMMON  /  TIM  /  TIMF1.DATE.D 

DIMFNSION  S(3,3) 

CALL  SPNATT(T,RASC,DECL) 

S(3,l)  =  COS(DICL^*COS(RASC) 

S(3,2)=C0S(DECL)*SIN(RASC) 

S(3,3)=SIN(DECL) 

CALL  SATVEC(T,PX,PY,PZ) 

PDOTS=PX*S (3,1) +PY*S ( 3 ,2  )+PZ*S (3 ,3 ) 

PSQ=PX**2+PY**2+PZ**2 

S1DEN=SQRT(PSQ~PD0TS**2) 

S(1,1)=(-PX+PD0TS*S(3,1) )/SlDEN 

S(1,2)=(-PY+PD0TS*S(3,2) )/SlDEN 

S(1,3)=(-PZ+PD0TS*S(3,3))/S1DEN 

S(2,l)=S(3,2)*S(l,3)-S(3,3)*S(l,2) 

S(2,2)=S(3,3)*S(1,1)-S(3,1)*S(1,3) 

S(2,3)=S(3,l)*S(l,2)-S(3,2)*S(l,l) 

RETURN 

END 

SIZE     145    00221 
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SUBROUTINE  SATVEC(T,PCX,PCY,PCZ) 
C      COMPUTES  THF  COMPONENTS  (PCX,PCY,AND  PCZ)  OF  THF  VECTOR  TO  THE 
C       SATELLITE  IN  FARTH  CENTERED  INERTIAL  COORDINATES  AT 
C         TIMF  (T)  BASED  ON  THE  CHEBYCHEV  COEFFICIENTS  CX(I ) , CY ( I ) ,CZ( I ) 

PIMFNSION  BX(13),BY(13),BZ(13) 

COMMON  /  SATPOS  /  CX ( 11 )  ,CY ( 11 ) , C Z ( 1 1 ) 

COMMON  /  TIM  /  TIME1,DATE,D 

CALL  NTIM(T,U) 

DO  5  1=12,13 

BX(I)=0.0 

BY(I)=0.0 

BZ(I)=0.0 
5   CONTINUE 

DO  15  J=l,ll 

BX(12-J)=CX(12-J)+2.0*U*BX(l2-J+l)-BX(12-J+2) 

BY(12-J  ^  =  CY(12-jU2.0*U*BY(12-J  +  l)-BY(12-J  +  2) 

BZ ( 12- J )=CZ ( 12-J ) +2 .0*U*BZ ( 12-J+l )-BZ( 12-J+2 ) 
15   CONTINUE 

PCX=(BX(1)-BX(3) )/2.0 

PCY=(BY(l)-BY(3))/2.0 

PCZ=(BZ(l)-BZ(3))/2.0 

RETURN 

END 

SIZE      156     00234 
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SUBROUTINF  SPNATT (T ,  SPRASC , SPDECL) 
C  COMPUTES  THE  RIGHT  AS  CENS  ION  (  SPRASC  )  AND  DECLINATION  (  SPDFCD 

C  OF  THE  SPIN  AXIS  AT  TIME  (T)  WITH  A  LINEAR 

C  INTERPOLATION  BASED  ON  VALUFS  OF  SPRA1.2  AND  SPDC1.2 

C 

C  THE  INTERPOLATION  TECHNIQUE  USED  HERE  IS  A  SIMPLE  APPROXIM- 

C  ATION  BUT  IS  SUFFICIENTLY"  PRECISE  FOR  THF  SMALL  PRECESSION 

C  RATES  ENCOUNTERED. 

C 

C  SPRA1,2  HAS  RANGE  -180  TO  +180  DEGREES 

C  SPDC1,2  HAS  RANGE  -90  TO  +90  DEGREES 

C  SPRASC  HAS  RANGE  -PI  TO  +PI  RADIANS 

C  SPDFCL  HAS  RANGE  -PI/2  TO  +PI/2  RADIANS 

C  NOTE   IF  THE  ABSOLUTE  VALUE  OF  THE  DIFFERENCE  BETWEEN 

C  SPRA1  AND  SPRA2  IS  GREATER  THAN  90  DEGREES,  THE 

C  INTERPOLATION  OF  DECLINATION  WILL  BE  INCORRECT 

COMMON/  SATATT  /  SPRA1 ,SPRA2,SPDC1 ,SPDC2 

COMMON  /  TIM  /  TIME1,DATE,D 

PI=3. 1415926535 

CALL  NTIM(T,U) 

RASC1=SPRA1*PI/180.0 

RASC2=SPRA2*PI/180.0 

DECL1=SPDC1*PI/180.0 

DFCL2=SPDC2*PI/180.0 

IF(RASC1.LT.0.0)  RASC1=2.0*PI+RASC1 

IF (RASC2.LT. 0.0)  RASC2=2.0*PI+RASC2 

SPRASC=0.5*(RASC2+RASC1+(RASC2-RASC1 )*U) 

SPRASC=AMOD(SPRASC,2.0*PI) 

IF(SPRASC .GT.PI )  SPRASC=SPRASC-2.0*PI 

SPDECL=0.5*(DECL2+DECL1+(DFCL2-DECL1)*U) 

RETURN 

END 

SIZE     101     00145 
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SUBROUTINE  NTIM(T.U) 
C  COMPUTFS  NORMALIZED  TIME  (U)  FROM  ACTUAL  TIME  ( T 1  BASED 

C  ON  TIME1  AND  D. 

C  TIME1  IS  IN  SECONDS 

C  D  IS  13  HOURS  IN  SFCONDS 

COMMON  /  TIM  /  TIME1,DATF,D 

IF  (TIMF1 .GT.T)  GO  TO  5 

U=2.0*(T-TIME1 )/D-l .0 

GO  TO  10 
5    U=2.0*(T  +  24. 0*60.  0*60.  0-TIMED/D-1.0 
10   RETURN 

END 

SIZE      33    00041 
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SUBROUTINE  GRANG(T.W) 

COMPUTE  GREENWICH  ANGLEU)  IN  RADIANS  FROM 
TIME(T)  BASED  ON  VALUES  OF  GRA1  AND  GRA2  ( 

RANGF  IS  -180  TO 
W  HAS  RANGE  -PI  TO  +PI  IN  RADIANS 

COMMON  /  GRNICH  /  GRA1 ,GRA2 
COMMON  /  TIM  /  TIME1,DATE,D 
PI=3. 1415926535 
CALL  NTIM(T,0) 
W1=GRA1*PI/180.0 
H/2=GRA2*PI/180.0 
IF(W1.LT.0.0)  W1=2.0*PI  +  Wl 
IFU2.LT. 0.0)  W2=2.0*PI  +  W2 
IF(W2.LT.tfl)  tf2=W2+2.0#PI 
W=0 . 5* ( W2+W 1+ ( W2-W1 ) *U ) 
W=AMOD(W,2.0*PI) 
IF(W.GT.PI)  W=W-2.0*PI 
RETURN 
END 


IN  DEGREES) 
+180  DEGREES 


SIZE 


84 
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NOAA  SCIENTIFIC  AND  TECHNICAL  PUBLICATIONS 

The  National  Oceanic  and  Atmospheric  Administration  was  established  as  part  of  the  Department  of 
Commerce  on  October  3,  1970.  The  mission  responsibilities  of  NOAA  are  to  assess  the  socioeconomic  impact 
of  natural  and  technological  changes  in  the  environment  and  to  monitor  and  predict  the  state  of  the  solid  Earth, 
the  oceans  and  their  living  resources,  the  atmosphere,  and  the  space  environment  of  the  Earth. 

The  major  components  of  NOAA  regularly  produce  various  types  of  scientific  and  technical  informa- 
tion in  the  following  kinds  of  publications: 


PROFESSIONAL  PAPERS  —  Important  definitive 
research  results,  major  techniques,  and  special  inves- 
tigations. 

CONTRACT  AND  GRANT  REPORTS  —  Reports 
prepared  by  contractors  or  grantees  under  NOAA 
sponsorship. 

ATLAS  —  Presentation  of  analyzed  data  generally 
in  the  form  of  maps  showing  distribution  of  rainfall, 
chemical  and  physical  conditions  of  oceans  and  at- 
mosphere, distribution  of  fishes  and  marine  mam- 
mals, ionospheric  conditions,  etc. 


TECHNICAL  SERVICE  PUBLICATIONS  —  Re- 
ports containing  data,  observations,  instructions,  etc. 
A  partial  listing  includes  data  serials;  prediction  and 
outlook  periodicals;  technical  manuals,  training  pa- 
pers, planning  reports,  and  information  serials;  and 
miscellaneous  technical  publications. 

TECHNICAL  REPORTS  —  Journal  quality  with 
extensive  details,  mathematical  developments,  or  data 
listings. 

TECHNICAL    MEMORANDUMS  —  Reports     of 

preliminary,  partial,  or  negative  research  or  technol- 
ogy results,  interim  instructions,  and  the  like. 


Information  on  avaitability  of  NOAA  publications  can  be  obtained  from: 

ENVIRONMENTAL  SCIENCE  INFORMATION  CENTER  (D822) 

ENVIRONMENTAL  DATA  AND  INFORMATION  SERVICE 

NATIONAL  OCEANIC  AND  ATMOSPHERIC  ADMINISTRATION 

U.S.  DEPARTMENT  OF  COMMERCE 

6009  Executive  Boulevard 
Rockville,  MD  20852 
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